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ABSTRACT 
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An  improved  version  of  a  previously  reported  long  bone  injury  criterion 
has  been  developed.  This  new  criterion  was  based  on  fitting  the  Ramberg- 
Osgood  equation,  through  optimization  techniques,  to  existing  stress-strain- 
strain  rate  curves  for  compact  bone  and  using  this  equation  to  solve  for  the 
strain  at  a  time  step  given  the  current  stress  and  the  strain  from  the  previous 
time.  Refinements  have  also  been  made  to  the  geometric  mqdeling  of  the  bones 
where  adequate  data  existed.  When  coupled  with  the^ATBj model  and  relations 
between  ultimate  stress  and  strain  rate,  this  criterion  has  provided  an 
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I .  INTRODUCTION 


The  Articulated  Total  Body  (ATB)  Model  has  been  developed  as  a  tool  to  predict 
and  analyze  the  response  of  the  human  body  to  potential  injury  causing  environments. 
This  model  has  been  primarily  used  for  the  analysis  of  vehicular  crashes,  and  has 
been  extreme ly  useful  in  predicting  gross  motion  of  the  body.  The  model  also  predicts 
the  loads  on  each  of  the  articulated  joints.  The  USAF  has  been  involved  in  the  devel¬ 
opment  and  testing  of  this  computer  model  and  has  applied  it  to  the  analysis  of  the 
response  of  pilots  to  ejection  from  jet  aircraft.  During  these  events  the  body  is 
subjected  to  high  accelerations  from  the  ejection  itself,  to  impacts  of  the  limbs 
with  each  other  and  with  the  cockpit,  and  to  "wind  flail"  from  the  high  velocity  wind 
stream  which  impacts  the  body  after  clearing  the  aircraft.  It  is  hoped  that  a 
better  understanding  of  these  events  will  lead  to  safer  designs  and  lower  injury 
rates. 

Despite  the  excellent  results  obtained  from  the  ATB  model,  the  current  version 
does  not  provide  the  information  necessary  to  predict  whether  the  event  will  or  will 
not  cause  an  injury.  It  is  therefore  necessary  to  augment  the  current  model  by  the 
addition  of  some  sort  of  injury  criterion  in  order  for  the  full  potential  of  the  ATB 
to  be  realized.  With  this  addition,  the  model  can  be  used  to  compare  different 
acceleration  profiles,  restraint  systems  and  other  variables  as  to  their  injury 
preventing  potential.  Currently  these  assessments  must  be  made,  in  a  very  qualitative 
way,  on  the  basis  of  motion  of  the  limbs  and  loads  in  the  joints. 

Since  the  majority  of  serious  injuries  resulting  from  ejections  are  bone 
fractures,  it  is  of  particular  interest  to  estimate  the  likelihood  of  long  bone 
fracture.  (It  should  be  noted  that  a  separate  computer  model,  the  Head-Spine  Model, 


is  specifically  designed  to  study  the  question  of  head  and  spine  injuries,  and  we 
have  consciously  ignored  this  aspect  of  injury). 

The  research  carried  out  under  this  minigrant  was  begun  by  the  principal 
investigator  during  a  Summer  Faculty  Research  Program  stay  at  WPAFB.  The  major 
objectives  of  that  research  program  were  to  establish  a  simplified  injury  criterion 
for  bone  and  to  implement  that  criterion  in  a  computer  program  (BREAK)  compatible 
with  the  existing  ATB.  Due  to  the  constraints  of  such  a  short  time  period  (10  weeks), 
the  project  was  necessarily  restricted  to  a  rather  narrow  scope.  Within  these  con¬ 
straints  the  objectives  were  accomplished.  A  detailed  report  of  this  project  can 
be  found  in  the  Final  Report  submitted  to  SCEEE,  and  published  as  a  technical  report 
by  AFAMRL  (1),  but  the  significant  developments  will  be  summarized  here. 

Injury  criteria  reported  in  the  literature  have  focused  on  two  distinct  approaches 
"experimental"  and  "analytical".  The  experimental  approach  comes  largely  from  the 
automobile  industry  and  involves  the  experimental  determination  of  loads  which, 
under  simulated  crash  tests,  cause  injury  to  cadavers.  For  long  bones  this  is 
restricted  to  the  Femur  Injury  Criterion  which  is  a  measure  of  the  axial  load 
applied  to  the  shaft  of  the  femur  which  causes  fracture.  The  current  Federal  Standard 
is  1700  lb  force  without  any  loading  rate  dependence.  A  number  of  authors  have 
suggested  different  criteria  with  loading  rate  sensitivity  (2,3).  Even  so,  the 
basic  fault  of  these  criteria  is  their  specificity  -  they  say  nothing  about  other  long 
bones  or  other  loading  conditions.  This  type  of  criterion  was  therefore  inappropriate 
to  a  study  of  ejection  seat  injuries. 

The  alternate  type  of  injury  criterion  is  based  on  the  analytic  approach  using 
the  material  and  geometric  properties  of  the  bone  and  the  calculated  stresses.  This 
approach  has  the  advantage  of  being  applicable  to  any  loading  situation  if  the 
properties  of  the  bone  are  known.  This  was  the  approach  taken. 
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The  material  properties  of  bone  are  highly  variable,  depend  on  a  large  number 
of  parameters  and  are  time  dependent  (4,5,6).  In  order  to  accurately  model  the 
behavior  of  bone,  consistant  data  including  complete  strain  rate  dependent  stress 
strain  curves  are  needed  for  human  bone.  The  data  which  most  nearly  fulfilled  these 
criteria  was  the  work  of  McElhaney  (7)  which  was  for  embalmed  human  bone  in  compres¬ 
sion.  His  results  indicated  the  following  relation  between  ultimate  stress  (a)  and 

O 

strain  rate  (e) 

a  =  4200  Log  e  +  33000.  (1) 

This  equation,  after  modifications,  became  the  basis  for  the  injury  criteria 
developed.  Two  major  changes  to  this  equation  were  made  -  the  strain  rate  dependency 
was  replaced  by  a  stress  rate  dependence,  and  the  constant  term  was  modified  to 
bring  the  results  in  line  with  published  static  results  for  fresh  human  compact  bone. 
The  rate  modification  was  accomplished  by  first  finding  a  relation  between  elastic 
modulus  and  strain  rate,  and  then  differentiating  the  elastic  stress  strain  relation 
with  respect  to  time  for  a  constant  strain  rate.  The  second  modification  was  based 

O 

on  e  =  .001  as  a  static  strain  rate  and  was  straightforward.  The  resulting  equation 
for  compression  (and  in  psi  units)  is 

a  =  3936  log  e  +  12000  (2) 

Similar  expressions  were  generated  for  tension  and  shear  based  on  the  assumption 
that  the  material  behavior  (rate  dependence)  would  be  the  same  in  the  other  modes 
and  that  only  the  constant  offset  would  change.  These  equations  formed  the  basis 
of  the  injury  criterion  models  developed  that  sumner. 

An  alternative  to  these  criterion  was  generated  using  stress  pulse  duration, 
rather  than  stress  rate,  as  the  time  variable.  This  alternate  approach  was  examined 
because  most  published  criteria  from  the  auto  industry  are  based  on  pulse  duration. 
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The  pulse  duration  was  approximated  by  assuming  a  sinusoidal  pulse  for  the  stress. 

The  "pulse"  criterion  tended  to  be  slightly  less  conservative  than  the  stress  rate 
criteria. 

The  stress  calculations  were  based  on  a  very  simplified  geometry,  namely 
straight,  uniform,  isotropic  hollow  tubes  for  all  the  long  bones.  Principal  stresses 
were  then  calculated  at  some  predetermined  number  of  sections  equally  spaced  along 
the  bone  axis  and  the  peak  tension,  compression  and  shear  stresses  were  monitored. 

The  peak  stresses  as  a  function  of  time  for  a  particular  limb  were  therefore  not 
necessarily  stresses  at  the  same  point,  the  position  along  the  axis  of  the  bone 
could  vary  with  time  also. 

A  test  case  was  run  based  on  an  ejection  simulation  to  demonstrate  the  program. 
Results  are  presented  in  the  Technical  Report  mentioned  above  and  a  sample  figure 
is  reproduced  here.  Figure  1  shows  the  maximum  compressive  stress  as  a  function 
of  time  for  the  left  lower  arm.  There  was  contact  at  .07  secs  and  some  high  level, 
higher  frequency  loading  at  the  elbow  near  the  end  of  the  data  used  (only  the  first 
200  msec  of  the  event  were  examined).  The  stress  rate  dependence  of  the  allowable 
stress  is  evident  as  is  the  slightly  less  conservative  nature  of  the  pulse  criterion. 
Note  that  the  pulse  criterion  has  only  been  calculated  for  the  maximum  stress  pulse 
and  hence  appears  as  a  constant. 

This  summarizes  the  state  of  this  project  at  the  start  of  the  minigrant. 
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COMPRESS.  STRESS  VS  ALLOWABLE 


II.  OBJECTIVES 


As  outlined  in  the  grant  proposal,  there  were  five  main  objectives  of  this 
investigation,  and  they  were: 

1.  Improve  the  efficiency  of  the  developed  program,  BREAK. 

2.  Search  out  and  reexamine  the  available  data  base  on  bone  material 
properties. 

3.  Develop  more  precise  constitutive  equations  for  bone. 

4.  Provide  a  statistical  analysis  of  the  available  data. 

5.  Redefine  the  injury  criteria  based  on  the  new  information. 


III.  RESULTS 


In  this  section  each  of  the  objectives  described  in  the  previous  section  will 
be  addressed  in  turn,  and  the  outcome  of  each  investigation  will  be  detailed. 

A.  The  greatest  problem  with  BREAK,  as  discussed  in  the  proposal,  has  been  the 
lack  of  compatibility  between  the  output  of  the  ATB  Model  and  the  input  needs  of 
BREAK.  Basically  the  problems  are  ones  of  transformation  of  coordinate  systems 
and  reconstruction  of  contact  forces  from  resultants.  As  an  example,  suppose  that 
a  limb  contacts  a  panel  resulting  in  normal  and  friction  forces.  The  output  of 
the  ATB  Model  will  include  the  linear  and  angular  position,  velocity  and  acceleration 
of  the  center  of  mass  of  the  segment  along  with  the  forces  and  moments  resulting 
at  the  joints  of  the  segment.  In  addition,  the  resultant  contact  force  and  the  point 
of  contact  will  be  given.  For  BREAK,  since  it  calculates  stress  levels  at  a  number 
of  points  along  the  axis  of  the  bone,  all  forces  must  be  known  in  local  coordinates 
and  the  contact  forces  must  be  known  in  a  vector  sense,  rather  than  in  magnitude 
only.  The  reconstruction  of  the  contact  forces  is  a  laborious  task  which  repeats 
the  work  already  done  by  the  ATB  program.  In  fact,  all  necessary  information  exists 
within  the  ATB  program  at  some  point  during  the  calculations  and,  as  proposed  earlier, 
a  file  could  be  set  up  containing  the  appropriate  information.  (See  program  CONTACT 
in  Appendix).  Unfortunately,  this  has  not  yet  been  carried  out.  The  major  difficulty 
has  been  the  lack  of  available  memory  on  the  CDC-CYBER  computer  system  used  by  AFAMRL. 
All  "non  essential"  data  are  currently  eliminated  during  processing  of  the  program  in 
order  to  save  space,  and  that  loss  currently  includes  the  data  which  is  needed  for 
BREAK.  Therefore  the  current  version  of  BREAK  (see  CONTACT  in  Appendix)  includes  all 
of  the  manipulations  which  must  be  performed  to  transform  the  given  data  into  the 
needed  format. 


Another  important  change  to  the  program  dealt  with  bone  geometry.  Rather  than 
continue  usinguniform,  hollow  tubes  to  represent  bones,  we  have  tried  to  imitate  some 


of  the  variations  in  actual  bone  geometry.  Using  the  limited  amount  of  available 
bone  geometry  data  (8,9)  we  have  been  able  to  model  variations  in  cross-section 
properties  along  the  bone  axes  (for  20  to  80%  of  distal  distance)  for  the  femur 
and  tibia.  No  data  was  found  for  the  upper  extremity  bones.  For  this  model  the  bone 
is  assumed  to  be  made  up  of  a  discontinuous  series  of-short  sections  of  uniform 
hollow  tubes. 

It  was  necessary  to  scale  the  available  data,  which  appeared  to  be  from  approx¬ 
imately  50th  percentile  subjects,  up  to  a  95th  percentile  man  so  that  the  bone  data 
would  coincide  with  the  other  pilot  data.  This  was  accomplished  by  making  the 
assumption  that,  for  different  size  bodies,  the  stress  levels  in  the  bones  would  be 
relatively  constant.  Using  the  body  weight  and  limb  length,  relative  ratios  of 
cross-sectional  area  and  area  moment  of  inertia  could  be  generated  using  axial 
compressive  stress  and  bending  stress  equations.  The  data  used,  in  both  original 
and  scaled  form,  are  shown  in  Figures  2  and  3.  See  program  STRESS  for  a  complete 
description  of  the  analysis,  along  with  the  calculation  of  stresses  at  various  cross- 
sections. 

B.  It  has  been  well  established  that  the  elastic  and  ultimate  properties  of 
compact  bone  are  loading  rate  dependent  (e.g.  7,10).  Ultimate  stress  levels  typically 
vary  by  a  factor  of  two  or  more  over  a  range  of  5  or  6  decades  of  strain  rates.  This 
effect  is  simply  too  large  to  ignore.  However,  the  use  of  ultimate  stress  values 
alone,  even  if  they  are  well  known,  is  Insufficient 

The  ATB  model  represents  the  body  as  a  series  of  rigid,  though  resilient,  links 
connected  by  springs.  As  such,  for  each  segment  the  end  loads  and  motions  are  well 
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FIGURE  3 


defined.  These  knowns,  along  with  the  bone  geometries,  are  sufficient  to  describe 
the  stress  state  at  any  timestep  and  at  any  point  along  the  limb.  The  strains  are 
then  related  to  the  stress  through  the  theory  of  elasticity  (for  an  elastic  structure). 
For  bone,  the  modulus  of  elasticity  is  a  function  of  strain  rate,  and  that  function 
is  not  analytically  invertable. 

Three  sets  of  stress-strain  curves  for  various  strain  rates  have  been  uncovered 

in  the  literature.  McElhaney  published  the  results  of  both  bovine  and  embalmed  human 

-3  -1 

bone  in  compression  using  constant  strain  rates  between  10  and  1500  sec  (7). 

More  recently.  Wood  has  reported  on  the  tensile  properties  of  fresh  cranial  bone  for 
strain  rates  from  .003  <  £  z  150  sec"^  (11).  This  data  is  of  little  utility  for  the 
ATB  model  because  Wood  used  very  small  bone  samples  and  because  cranial  bone  is 
significantly  different  in  structure  from  diaphysial  compact  bone.  Crowninshield 
and  Pope  (12)  have  also  published  results  for  compact  bovine  bone  in  tension  over 
strain  rates  of  .00167  s  t  <  250.  Their  results  indicate  a  much  larger  plastic  strain 
region  for  longitudinal  samples  than  do  the  results  previously  reported.  However, 
Burstein  et.al_.  (13)  have  also  shown  considerable  plastic  strain  in  bovine  bone  in 
tension  and  very  little  in  compression.  In  Burstein's  tests  strain  rates  were  not 
reported,  but  loading  duration  was  between  1/10  and  1/2  sec. 

The  above  three  papers  constitute  a  very  small  data  set  and  none  of  the  results 
are  for  fresh  human  long  bone.  Other  researchers  have  dealt  with  portions  of  this 
problem,  however  they  usually  report  only  ultimate  properties  and  not  the  full  stress- 
strain  curves.  Lewis  and  Goldsmith  (14),  for  example,  used  a  split  Hopkinson  Bar 
method  to  measure  the  fracture  stress  of  bovine  bone  in  compression.  The  strain 
histories  were  pulses  rather  than  constant  strain  rate.  The  fracture  strain  rates 
were  calculated  by  dividing  the  fracture  strain  by  the  time  and  these  rates  varied 
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over  approximately  6  orders  of  magnitude.  These  results  are  consistently  higher 

than  the  results  of  McElhaney.  Wright  and  Hayes  (15)  reported  ultimate  tensile 

-4  •  -1 

strength  for  fresh  bovine  bone  for  5.3  x  10  i  e  <  237  sec  .  They  also  present 
two  characteristic  load-displacement  curves,  but  no  stress-strain  data. 

C.  Much  progress  has  been  made  in  the  area  of  establishing  analytic  functions 
which  describe  the  relationship  between  stress,  strain  and  strain  rate.  Using 
available  published  data  and  our  optimizing  curve  fit  computer  program  we  have  been 
able  to  get  good  fits  of  the  data  using  a  standard  Ramberg -Osgood  equation, 

e  =  +  aoN£b  .  (3) 

ce 

This  equation  is  a  standard  equation  for  modeling  visco  elastic-plastic  behavior 
of  materials  (16).  Results  for  the  three  available  sets  of  curves  are  shown 
numerically  in  Table  1,  and  graphically  in  Figures  4,  5,  and  6.  The  ability  of 
this  function  to  fit  such  a  wide  variation  of  results  is  encouraging. 

Because  of  the  orders  of  magnitude  variations  in  the  strain  rate,  certain 
constants  in  equation  (3)  were  extremely  sensitive.  For  example,  the  constant 
multiplier  "a"  for  the  plastic  term  varied  by  as  much  as  55  orders  of  magnitude 
from  author  to  author.  It  was  therefore  necessary  to  take  our  optimization  approach 
(see  program  MYFIT  in  Appendix)  to  minimize  the  function  F  (the  sum  of  weighted 
errors  of  data-point  fits),  where  F  is  defined  by  equation  (4). 

M  M 

F  =  E  [H(e  -  e0)]N  +  r  E  <gi>z  (4) 

i=l  p  0  i=l 


12 


STRESS  -  STRAIN  CURVES 

BOVINE  :  TENSION  (CROWNINSHlELD  &  POPE) 
LONGITUDINAL 


T~ T ■■■■■■■■  r  1  r~p-  'I  ■  t 

006  0.012  0.018  0.024  0.030  0.036  0.042  0.048 

STRAIN  (IN/IN) 


FIGURE  5 


STRESS  -  STRAIN  CURVES 

CRANIAL  COMPACT  80NE  :  TENSION  (WOOD) 


0.0030  0 . 0045 


0.001S 


STRAIN  (IN/ IN) 

FIGURE  6 


0.0060 


and  where 

M  =  the  number  of  experimentally  observed  points  of  (a,e,e) 

£p  =  strain  predicted  by  curve 
eQ  =  strain  observed  by  experiment 
r,H  =  constant  multipliers  for  weighting 
N,z  =  constant  powers  for  weighting 
<gi>  =  exterior  penalty  function  for  weighting 

The  objective  function  F  is  minimized  in  each  of  the  five  coordinate  directions 
(a,b,N,d.c)  using  a  combination  of  golden  section  search  and  parabolic  interpolation 
methods  (17).  The  resultant  coordinates  are  the  coefficients  of  the  equation 
that  best  fits  the  set  of  observed  data. 

The  results  of  the  optimization  process  are  the  curves  shown  in  Figures  4,  5 
and  6.  It  can  be  seen  that  there  is  excellent  agreement  between  most  of  the  predicted 
and  experimental  results.  The  largest  variations  occur  in  the  very  high  strain  rate 
tests  of  McElhaney.  Forcing  the  function  to  try  and  fit  the  1500  sec”^  data  tends 
to  also  skew  the  rest  of  the  curves.  As  an  experiment,  the  1500  sec-1  data  was 
removed  from  the  sample  population  and  the  curves  refit  to  the  remaining  data.  This 
trial  showed  a  very  good  correlation  with  the  remaining  data  and  a  high  prediction 
at  1500  sec"^ .  As  this  is  the  only  data  available  at  this  high  a  strain  rate  it 
is  difficult  to  draw  any  conclusions  about  this  response. 

D.  It  is  obvious  from  the  previous  section  that  no  meaningful  statistical 
analysis  can  be  performed  on  the  meager  amount  of  available  data.  However,  it  is 
interesting  to  compare  the  various  results  and  this  is  shown  in  Figures  7,  8,  9 
and  10.  These  figures  show  predicted  stress-strain  curves  at  a  given  strain  rate 
based  on  various  authors'  results.  Be  aware  first  of  all  that  these  curves,  except 
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STRESS  -  STRAIN  CURVES 

STRAIN  RATE  =  100  /  SEC 


STRAIN  (IN/IN) 


STRESS  -  STRAIN  CURVES 

STRAIN  RATE  =  1  /  SEC 


17 


0,00  0.01  0.02 


STRESS  -  STRAIN  CURVES 

STRAIN  RATE  =  .01  /  SEC 


STRESS  -  STRAIN  CURVES 

STRAIN  RATE  =  .001  /  SEC 


for  the  two  by  Wood,  all  represent  different  tests -embalmed  human  femur  in  compres¬ 
sion;  human  cranial  compact  bone  in  tension;  and  bovine  femur  in  tension  for  long¬ 
itudinal  and  transverse  directions.  In  that  sense  we  are  comparing  apples  to  oranges. 
For  this  reason,  these  figures  are  shown  only  for  general  comparison.  It  had  been 
hoped  that  a  more  direct  comparison  of  similar  data  from  the  literature  could  have 
been  made. 

The  cranial  bone  samples  are  generally  weaker  with  stiffnesses  bracketed  by 
the  compression  results  (stiffest  and  highest  strength)  and  the  tension  results  (less 
stiff  and  most  plastic).  Results  for  transverse  tension  samples  are  shown  to  be 
significantly  lower  in  strength  and  stiffness.  As  mentioned  earlier,  the  only 
significant  plastic  zone  was  found  by  Crowninshield  and  Pope  for  bovine  bone  in 
longitudinal  tension.  The  results  of  Wood  indicate  that  there  is  some  consistency 
among  data  from  tests  on  similar  material  and  under  similar  conditions. 

A  study  was  carried  out  during  the  initial  summer's  work  to  assemble  a  small 
subset  of  this  data,  namely  the  ultimate  strengths  for  fresh  human  compact  bone  from 
femurs  tested  statically  in  the  longitudinal  direction.  The  results  are  reproduced 
in  Table  2.  For  the  four  modes  of  failure,  the  means  and  standard  deviations  are, 
respectively;  tension  99.3/26/3,  compression  174.4/34.2,  shear  69.0/14.8  and  bending 
161.4/11.9,  all  *10**6  N/m**2.  These  figures  illustrate  the  wide  discrepencies 
which  exist  for  even  a  narrowly  defined  subset  of  available  data.  The  wide  variations 
in  these  values  could  be  due  to  a  number  of  causes  including  the  age  and  moisture 
content  of  the  specimens  and  the  test  procedures  utilized  by  the  investigators. 

E.  The  injury  criteria,  as  first  developed,  were  based  on  relationships  between 
the  ultimate  stress  and  the  strain  rate  (equation  1).  This  equation  was  then  modified, 
based  on  a  constant  strain  rate  approximation,  to  yield  a  relation  between  ultimate 


Table  2 


Ultimate  Strengths  of  Fresh  Human  Compact  Bone  From  Femurs,  Tested 
Statically  and  Stressed  in  the  Longitudinal  Direction 

{xlO6  N/m2/xl03  psi) 


Tension 

Compression 

122/17.7* 

159/23.1* 

86.5/12.5* 

193/28.0* 

133/19.3* 

134.5/19.5^ 

76. 2/11. 0+ 

210. 9/30. 6+ 

78. 9/11. 4+ 

Shear 

Bendinq 

53.1/7.7* 

152/22.0& 

82.4/11.9* 

153/22. 1* 

71.6/10.4* 

157/22.8* 

164/23.8* 

181/26.2* 

Compiled  from  various  sources  reported  in  (*)  Reilly  and  Burstein  (6),  (+)  Evans 
(4),  ($)  Reilly  and  Burstein  (18),  (&}  Mather  (19)  and  (<£)  Vose  and  Kubala  (20). 
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stress  aid  stress  rate  which  could  be  calculated  directly  from  the  available 
stress  time  histories  (equation  2).  This  solution  is  rather  simplistic  and  not 
entirely  satisfactory.  A  more  systematic  approach  based  on  the  Ramberg -Osgood 
equation  (equation  3)  has  yielded  more  promising  results. 

The  elemental  question  is,  how  to  calculate  strain  when  only  the  stress  time 
history  is  known  and  the  relation  between  stress  and  strain  is  strain  rate  dependent? 

It  is  clear  on  inspection  that  equation  (3)  does  not  invert  or  separate  into 
convenient  parts,  particularly  since  the  factors  b  and  d  are,  in  general,  not  integers. 
The  approach  taken  has  been  one  of  a  numerical  approximation  based  on  the  following 
equation  at  time  step  i: 


£i  = 


c  e<-£i-l 


(^ 


M  ^4 ^ 
_i_  ,  _  N/  1  1  “  1  \ 

"d  aai  (  W~> 


At 


(5) 


Given  the  stress  at  a  time  step  and  the  strain  at  the  previous  time  step, 
equation  (5)  can  be  solved  numerically  for  the  current  strain.  This  solution  procedure 
is  easily  initialized  by  assuming  zero  stress  and  strain  at  time  =  0.  Solution  of 
this  equation  has  been  successfully  accomplished  and  Figure  11  shows  an  example 
stress  strain  response  to  the  sinusoidal  stress  wave  shwon  in  Figure  12.  The  solution 
displays  the  expected  hysteresis  type  response. 

There  were  a  number  of  stumbling  blocks  which  had  to  be  overcome  in  order  to 
achieve  this  solution.  The  second  term  of  equation  (5),  which  accounts  for  the  plastic 
portion  of  the  response,  can  have  a  severe  dest*'-'’!  < zing  effect  on  the  solution 
procedure,  depending  on  the  values  of  the  constants  and  the  strain  rate.  It  was 
therefore  necessary  to  use  a  two  step  method  based  on  the  Regula-Falsi  technique. 

An  initial  solution  was  generated  using  only  the  first  term  of  the  equation  to  be 
used  as  a  starting  point.  With  this  initial  guess  as  a  basis  the  final  solution 


22 


was  obtained  using  both  terms.  This  two  step  procedure  greatly  reduced  the  incidence 
of  non  convergence,  but,  even  using  double  precision,  did  not  entirely  eliminate 
these  difficulties  if  the  stress  levels  were  too  high. 

A  second  problem  was  encountered  due  to  the  relaxation  phase  of  the  loading. 

Since  b  and  d  are  not  integers,  a  negative  strain  rate  raised  to  such  a  power  has 
no  meaning.  Furthermore,  it  can  be  assumed  from  standard  hysteresis  loop  response  that 
relaxation  is  essentially  an  elastic  phenomenon.  For  these  reasons,  only  the  first 
term  of  the  equation,  with  appropriate  absolute  values,  is  used  during  relaxation. 

The  final  problem  arose  from  the  plastic  deformation  which  may  occur.  Since 
equation  (5)  presumes  a  zero  strain  at  zero  stress,  an  appropriate  plastic  strain 
offset  must  be  added  to  the  solution  after  any  relaxation.  This  plastic  effect  is 
clearly  shown  in  Figure  9. 

Using  the  solutions  generated  by  equation  (5)  it  is  now  possible  to  follow  the 
strain  and  strain  rate  response  of  the  system  and  to  directly  use  equation  (1) 
to  predict  ultimate  stress.  This  procedure  enables  one  to  realistically  use  available 
ultimate  stress  values  as  a  fracture  criterion. 

The  point  may  be  made  that  yield  rather  than  ultimate  stress  would  be  a  more 
conservative,  and  more  meaningful  injury  criterion.  Yield  stress  would  unquestionably 
be  more  conservative,  though  some  questions  must  first  be  answered.  First,  of  course, 
yield  point  must  be  defined.  This  definition  is  arbitrary,  depending  on  the 
material.  For  most  engineering  metals,  yield  is  defined  in  terms  of  that  stress  which 
produces  a  plastic  strain  of  .002.  There  is  no  such  accepted  standard  for  bone  and 
all  figures  reported  in  the  literature  are  for  ultimate  strengths.  If  a  yield 
point  definition  could  be  agreed  upon  for  bone,  then  yield  strengths  could  be  defined 
from  the  solution  of  equation  (5). 
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A  perhaps  more  interesting  question  is,  what  is  the  biological  effect  in  living 
bone  of  exceeding  the  "yield  strength?"  No  known  investigation  has  been  made  into 


this  phenomenon  and  its  relation  to  injury.  In  fact,  the  overall  level  of  under¬ 
standing  of  bone  fracture  mechanics  is  quite  low.  For  example,  it  is  conceivable 
that  reconstruction  in  living  bone  would  overcome  yielding  effects  without  injury. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


1 


Given  the  available  bone  fracture  data  and  the  constraints  of  the  current 
ATB  model,  a  reasonable  mechanism  has  been  developed  for  utilizing  a  long  bone 
fracture  injury  criterion  for  the  ATB  model.  Specifically,  this  mechanism  is 
based  on  a  mathematical  representation  of  the  behavior  of  bone  using  the  Ramberg- 
Osgood  equation  and  on  an  approximation  of  the  stress-time  history  of  the  bone 
generated  using  the  ATB  joint  and  contact  forces.  The  specific  injury  criterion 
which  should  be  employed  is  still  under  consideration.  Options  include:  ultimate 
strength,  which  would  be  non  conservative  but  for  which  considerable  data  exist; 
yield  strength,  which  would  be  conservative,  but  for  which  a  definition,  and  hence 
data,  do  not  exist;  some  fraction  of  either  of  these  as  a  factor  of  safety;  or  some 
other  criterion  such  as  maximum  strain  or  strain  energy.  At  present  the  ultimate 
strength,  reduced  by  a  factor  of  safety,  is  thought  to  be  the  best  criterion.  (The 
pulse  criterion  discussed  in  the  first  report  (1)  is  felt  to  have  value  only  as  a 
comparison  to  existing  automobile  standards  and  hence  to  have  no  applicability  to  the 
pilot  ejection  problem.)  The  resulting  injury  criterion  program  remains  in  the 
"post  processor"  mode,  i.e.  it  is  not  integrated  into  the  ATB  model,  and  this  is  the 
preferred  configuration  until  the  methodology  becomes  more  fixed. 

The  Ramberg -Osgood  equation  appears  to  provide  a  very  good  means  of  mathematic¬ 
ally  modelling  the  stress-strain-strain  rate  behavior  of  bone.  The  fact  that  this 
equation,  through  proper  optimization  techniques,  can  be  used  to  represent  the  three 
very  different  responses  of  bone  shown  here  from  the  literature  leads  one  to  suspect 
that,  whatever  the  "true"  response  of  bone  may  be,  it  can  be  adequately  described  by 
equation  (3).  This  flexibility,  coupled  with  the  dearth  of  good  data  on  dynamic 
bone  properties,  was  a  prime  motivating  factor  in  taking  this  approach. 
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Further  areas  of  inquiry  point  strongly  in  the  direction  of  establishing  a 
more  comprehensive  and  consistent  data  base.  The  state-of-the-art  in  modeling 
has  clearly  outstripped  both  the  available  information  on  bone  material  properties 
and  whole  bone  geometry,  and  the  understanding  of  the  mechanisms  of  bone  fracture. 

A  thorough  understanding  of  these  properties  and  mechanisms  is  necessary  before  any 
strides  can  be  made  toward  formulating  a  more  effective  long  bone  injury  criterion. 
The  present  criterion  is  designed  to  adapt  to  new  information  as  it  becomes  available 
through  modification  of  the  constants  in  the  Ramberg -Osgood  equation  and  represents, 
in  its  present  form,  the  most  sophisticated  injury  criterion  available. 
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7  :  OUTPUTS  LIST  OF  OBSERVED  AND  PREDICTED  VALUES  AT  TERMINATION. 
I/O  REQUIREMENTS  J 


*00000510 

*00000520 

*00000530 

*00000540 

*00000550 

*00000560 

*00000570 

*00000580 

*00000590 

*00000600 

*00000610 

*00000620 

*00000630 

*00000640 

*00000650 

*00000660 

*00000670 

*00000680 

*00000690 

*00000700 

*00000710 

*00000720 

*00000730 


*00000740 

FILE  *1  :  CONTAINS  INPUT  VARIABLES  AS  DESCRIBED  ABOVE  FOLLOWING  SDATA* 00000750 
CARD  AT  END  OF  PROGRAM.  *00000760 

FILE  #4  :  CONTAINS  OBSERVED  VALUES.  THE  FIRST  RECORD  MUST  BE  A  TITLE*00000770 
FILE  #5  :  CONTAINS  X  S  COEFF  ON  TWO  RECORDS.  UPDATED  AT  TERMINATION. *00000780 

*00000790 

**********************************■*<****************************♦*******00000800 


S OPTIONS  CCOHP=0 

ON  ERROR  GO  TO  70 
IMPLICIT  REAL  *  8  (A-H,  O-Z) 

REAL  NSCALE  (5) 

DIMENSION  X  (5)  ,  STRESS  (60)  ,  Y(60),  RATE  (60),  YAPROX  (60) 
DIMENSION  A  (5)  ,  DIFF  (60) 

INTEGER  TITLE  (20) 

COMMON  STPSSS,  Y,  RATE,  YAPROX,  M,  SSQ,  DIFF,  TIMES 
COMMON  /WEIGHT/  R,  Z,  ERROR,  NUHPEN,  PENAL 
COMMON  /SEARCH/  EPS 

COMMON  /FUN  1/  COEFF(5),  TOTAL,  UP,  A,  I 
EXTERNAL  FUNCT 

DATA  NSCALE  /0.90,  0.90,  1.00,' 1.00,  1.00  / 


SET  PROGRAM  PARAMETERS 


Z  »  2.0D0 
TOL  »  i.OD-7 
TIMES  »  1.0D3 


UP  =  2.0D0 
ERROR  =  1.0D-4 


R  *  10.0D0 
HAXFN  —  20 


00000810 

00000820 

00000830 

00000840 

00000850 

00000860 

00000870 

00000880 

00000890 

00000900 

00000910 

00000920 

00000930 

00000940 

00000950 

00000960 

OOOOC-70 

00000990 

00000990 

00001000 


uou 


CALCULATE  SQUARE  ROOT  OF  MACHINE  EPSILON  FOR  SUBROUTINE  FMIN 

EPS  =  1.0D00 
10  EPS  =  EPS  /  2.0DO 
TOL 1  =  1.0D0  ♦  EPS 
IF  (TOL 1  .GT.  1.0D00)  GO  TO  10  ' 

EPS  =  DSQRT  (EPS) 

KOUNT  =  NUHPEN  =  0 
READ  (1,*)  N 

READ  (5,*)  (1(1)  ,  I  »  1,  N) 

READ  (5,*)  (COEFF(I),  I  =  1,  N) 

C5  HEAD  (1,*)  (X  (I) ,  I  =  1,  N) 

C5  -BEAD  (1,*)  (COEFF(I),  I  *  1,  N) 

READ  (4,95)  TITLE 
S  =  1 

20  READ  (4,*,END=25)  RATE(H)  ,  Y  (H)  ,  STRESS<») 

a  =  a  i 

GO  TO  20 
25  a  =  a  -  i 

1*1;  POINT  =  X(I) 

WRITE  (3,98)  TITLE 

WRITE  (3,105)  TOL,  UP,  R,  2,  ERBOB 

WRITE  (3,106)  TIMES 

WRITE  (3,108)  (NSCALE(J),  J=  1,  N) 

WHITE  (3,107)  B 
DO  30  J  *  1,  N 

WRITE  (3,110)  J,  X(J),  COEFF(J) 

30  A  (J)  =  X(J)  *  COEFF(J) 

HOLD  =  FONCT  (POINT) 

WRITE  (3,100)  HOLD 

IF  (NUHPEN  .GT.  0)  WRITE  (3,125)  NUHPEN,  PENAL 
DO  60  K  *  1,  HAXFN 
IOUT  *  KOUNT 2  =  0 
DO  50  I  =  1,  N 
DO  40  J  »  1,  N 

40  A ( J)  =  X(J)  *  COEFF(J) 

C4  WRITE  (3,110)  I,  X(I),  COEFF(I) 

IF  (X(I)  .EQ.  0.0D0  .OR.  NSCALE(I)  .EQ.  0-0)  THEN  DO 
G  *  X(I) 

KOUNT2  *  KOUNT2  +1 
GO  TO  50 
ENDIF 

IF  (X(I)  -GT.  0.0D0)  THEN  DO 
GLOW  =  X(I)  *  (1.0  -  NSCALE  (I)  ) 

GHIGH  *  X(I)  *  (1.0  ♦  NSCALE  (I)) 

IF  (GLOW  .LT.  1.0D-3  .AND.  GLOW  . GT.  0.0D0)  GLOW  =  -GLOW 
ELSE  DO 

GLOW  =  X(I)  *  (1.0  ♦  NSCALE  (I)  ) 

GHIGH  *  X(I)  *  (1.0-  NSCALE  (I) ) 


00001010 

00001020 

00001030 

00001040 

00001050 

00001060 

00001070 

00001080 

00001090 

00001100 

00001110 

00001120 

00001130 

00001140 

00001150 

00001160 

00001170 

00001180 

00001190 

00001200 

00001210 

00001220 

00001230 

00001240 

00001250 

00001260 

00001270 

00001280 

00001290 

00001300 

00001310 

00001320 

00001330 

00001340 

00001350 

00001360 

00001370 

00001380 

00001390 

00001400 

00001410 

00001420 

00001430 

00001440 

00001450 

00001460 

00001470 

00001480 

00001490 

00001500 


IF  (GHIGH  . GT.  -1.0D-3  .AND.  GHIGH  . 

LT.  0.0D0) 

GHIGH  =  -GHIGH 00 00 151 0 

ENDIF 

00001520 

C4 

WRITE  (3,115)  GLOW,  GHIGH 

00001530 

G  =  FMIN  (GLOW,  GHIGH,  FUNCT,  TOL) 

0  001540 

TOTAL  =  FUNCT  (G) 

00001550 

C4 

WRITE  (3,120)  G,  SSQ 

00001560 

C4 

IF  (NUHPEN  .GT.  0)  WRITE  (3,125) 

NOMPEN,  PENAL 

00001570 

C4 

IF  (TOTAL  .GT.  HOLD)  THEN  DO 

00001580 

C4 

PRINT, 'BAD  STEP?' 

00001590 

C4 

ENDIF 

00001600 

C  4 

WRITE  (3,160)  TOTAL 

00001610 

45 

IF  (DABS  (1(1)  -G)  .GT.  TOL)  IOUT  = 

1 

00001620 

IF  (DABS  (  (TOTAL-HOLD) /HOLD)  .LT. 

TOL) 

KOUNT2  = 

KOUNT 2  +  1  00001630 

HOLD  =  TOTAL 

00001640 

50 

1(1)  =  G 

00001650 

C 

00001660 

C 

NORMALIZE  X'S  TO  BETWEEN  1.0  S  10.0.  CORRECT  CHANGE 

IN  VALUE  OF  00001670 

C 

COEFF  SO  PRODUCT  IS  SAME. 

00001680 

c 

00001690 

DO  55  I  =  1,  N 

00001700 

IF  (X(I)  .EQ.  0.0D0)  GO  TO  55 

00001710 

WHILE  (DABS  (X  (I))  .LT.  1.0D0)  DO 

00001720 

X(I)  =  X(I)  *  10.0D0 

00001730 

COEFF  (I)  =  COEFF  (I)  /  10.0D0 

00001740 

END  WHILE 

00001750 

55 

CONTINUE 

00001760 

c 

00001770 

c 

CHECK  AND  UPDATE  CONVERGENCE  CRITERIA  - 

TERMINATE  IF 

MET  00001780 

c 

00001790 

IF  (IOUT  .EQ.  0)  KOUNT  =  KOUNT  +  1 

00001800 

IF  (KOUNT.  EQ.  2)  THEN  DO 

00001810 

WRITE  (3,135)  TOL 

00001820 

WRITE  (3,145)  K 

00001830 

GO  TO  70 

00001840 

ENDIF 

00001850 

IF  (KOUNT2  .EQ.  N)  THEN  DO 

00001860 

WRITE  (3,140)  TOL 

00001870 

WRITE  (3,145)  K 

00001880 

GO  TO  70 

00001890 

ENDIF 

00001900 

c 

00001910 

60 

CONTINUE 

00001920 

WRITE  (3,147)  BAXFN 

00001930 

70 

DO  80  J  *  1,  N 

00001940 

80 

WRITE  (3,110)  J,  X(J),  COEFF  (J) 

00001950 

REWIND  5 

00001960 

WRITE  (5,130)  ,  (X(I)  ,  I  *  1,  N) 

00001970 

WRITE  (5,130)  (COEFF  (I),  I  *  1,  H) 

00001980 

WRITE  (3,150)  TOTAL 

00001990 

IF  (NUHPEN  .GT.  0)  WRITE  (3,155)  NUHPEN 

00002000 

i  B3  =  PENLTY  (1) 

'  WRITE  (3,165) 

r  DO  85  I  =  1,  (I 

r  WRITE  (3,170)  I,  STRESS  (I)  ,  RATE  (I)  ,  Y  (I)  ,  YAPROX  (I) 
f  05  CONTINUE 
STOP 

95  FOR 3 AT  (20A4) 

98  FORMAT  (1  HI, ‘TITLE  :  »,20A4/) 

100  FORMAT  ( 1 H  ,  *  AT  START,  FUNCT  =',  1PD21.13//) 

105  FORMAT  (1H-, 'PROGRAM  PARAMETERS T24, *  TOL  =' , 1PD9.2,/,T4, 

*  'ABS  ERROR  RAISED  TO  POWER: », 1PD9. 2,//, •  PENALTY  CONTROLS: »/,T 15 

*  'MULTIPLIER  =', 1PD12.5,/,T20, 'POWER  =• , 1PD9. 2 ,/,T14 , 

*  'ERROR  LEVEL  =' , 1PD12. 5//) 

106 •  FORMAT  (1H  , 'TIMES  =',1PD10-3) 

107  FORMAT  (1H  ,13,'  POINTS  WERE  OBSERVED'//) 

108  FORMAT  (*  SCALE  VECTOR  IS:  ',6P6.2/) 

110  FORMAT  {•  - X ( ' ,12, ')  *»  , 1PD21 . 13, •  *  *,1PD8.1) 

115  FORMAT  (»  SEARCH  RANGE  »,1PD12.5,'  TC  *,1PD12.5) 

120  FORMAT  (1H  , 'MINIMUM  POINT  3  G  =  ' , 1PD21 . 13,/,T 15, 

*  ' SSQ  =',  1PD21. 13) 

125  FORMAT  (I4,  «  CONSTRAINTS  VIOLATED  —  PENALTY  =',1PD13. 5) 

130  FORMAT  (5D21. 13) 

135  FORMAT  ('  CONVERGENCE  OF  ALL  X* 'S  TO  WITHIN* ,  1PD12. 3) 

140  FORMAT  ('  DELTA  FUNCT  HAS  BEEN  LESS  THAN* , 1PD 12.5, •  FOR  5  STEPS') 
145  FORMAT  (1H-, 'FINAL  RESULT  REACHED  AFTER  ',13,'  SERIES') 

147  FORMAT  flH-,'NO  CONVERGENCE  AFTER  *,I3,'  SERIES') 

150  FORMAT  (1H-,»AT  FINISH,  FUNCT  =  ',1PD21.13) 

155  FORMAT  ('OAT  FINISH,  ',13, »  CONSTRAINTS  VIOLATED:'/, 

*  1 H  ,T6,'#»  , T 10, 'STRESS  (KSI)  '  ,T28, 'RATE*  ,  ^ 

*  T38 , ' PRED.  STRAIN* ,T52,*OBS.  STRAIN* ,T69, 'DIFF'/) 

160  FORMAT  ('-FUNCT  =• , 1PD21. 1 3/) 

165  FORMAT  (1H0,/,T4, '#  *,T15, 'STRESS  (KSI) • ,T34, 'RATE* ,T45, 'OBS.  STRAIN 

*  ,T60, 'PRED.  STRAIN'/) 

170  FORMAT  (1H  ,T3,I2,T10, 4(1PD15.3) ) 

END 

DOUBLE  PRECISION  FUNCTION  FUNCT  (POINT) 

IMPLICIT  PEAL  *  8  (A-H,  O-Z) 

DIMENSION  X  (5)  ,  STRESS  (60)  ,1(60)  ,RATE(60)  ,YAPROX  (60)  ,DIFF(60) 
COMMON  STRESS,  Y,  RATE,  YAPROX,  M,  SSQ,  DIFF,  TIMES 
COMMON  /FUN  1/  COEFF(5),  TOTAL,  UP,  X,  II 
SSQ  =  0.0D0 

X(II)  *  POINT  *  COEFF(H) 

CALCULATE  PREDICTED  STRAINS  AND  SUM  OF  ERRORS. 

DO  10  I  *  1,  M 

IF  (STRESS (I)  .EQ.  0.0D0)  THEN  DO 
YAPROX (I)  *  DIFF (I)  *  O.ODC 
GO  TO  10 
RNDIF 


STEPS') 


00002019 

00002020 

00002030 

00002040 

00002050 

00002060 

00002070 

00002080 

00002090 

00002100 

,00002110 

00002120 

00002130 

00002140 

00002150 

00002160 

00002170 

00002180 

00002190 

00002200 

00002210 

00002220 

00002230 

00002240 

00002250 

00002260 

00002270 

00002280 

00002290 

00002300 

00002310 

•00002320 

00002330 

00002340 

00002350 

00002360 

00002370 

00002380 

00002390 

00002400 

00002410 

00002420 

00002430 

00002440 

00002450 

00002460 

00002470 

00002480 

00002490 

00002500 


TEMP2  »  STRESS  (I)  /  X  (5)  /  RATE  (I)  **  X(4) 

T'lMpa  =  RATE  (I)  **  X  (2)  *  STRESS  (I)  **  X{3) 

YAPROX  (I)  =  TEMP2  ♦  X  ( 1)  *  TEMP3 

APPLY  WEIGHTING  FACTORS  TO  RESIDUAL 

DIF  =  DABS (YAPROX (I)  -  Y(I)) 

SSQ  =  SSQ  ♦  (DIF  *  TIMES)  **  UP 
DIFF(I)  =  DIF 
10  CONTINUE 

TOTAL  =  SSQ  +  PSNLTY  (0) 

FUNCT  =  TOTAL 

3  -WRITE  (3,100)  POINT,  FUNCT 
RETURN 

100  FORMAT  ( 1 H  ,*AT  FUNCT,  POINT  =' , 1PD21. 13,  *  FUNCT  =* , 1PD2 1. 13) 
END 

DOUBLE  PRECISION  FUNCTION  PENLTY  (IDUM) 

IMPLICIT  REAL  *  8  (A-H,  0-2) 

DIMENSION  STRESS  (60)  ,Y{60)  , BATE  (60)  ,YAPROX(60)  ,DIFF(60) 

COMMON  STRESS,  Y,  RATE,  YAPROX,  M,  SSQ,  DIFF,  TIMES 
COMMON  /WEIGHT/  R,  2,  ERROR,  KOUNT ,  HOLD 
HOLD  =  0.0D0 
KOUNT  *  0 

CALCULATE  PENALTY  IF  ABS.  ERROR  GREATER  THAN  SPECIFIED  MAX. 

DO  101  -  1 ,  M  \ 

B3  =  DIFF (I) 

IF  (B3  .GT.  ERROR)  THEN  DO 
HOLD  =  HOLD  ♦  (B3  *  TIMES)  ♦»  2 
KOUNT  *  KOUNT  ♦  1 
IF  (IDUM  .EQ.  1)  THEN  DO 

WRITE  (3,100)  KOUNT, STRESS  (I)  , RATE  (I)  ,  YAPROX  (I)  ,Y(I)  ,B3 
ENDIP 
END  IF 

10  CONTINUE 

PENLTY  *  HOLD  =  HOLD  *  R 
100  FORMAT  ( 1 H  ,15,5  ( 1PD1 4.  5)  ) 

RETURN  ;  END 
SOPTIONS  NOLIST 

DOUBLE  PRECISION  FUNCTION  FHIN ( AX, BX, F,TOL) 

DOUBLE  PRECISION  AX,BX,?,T0L 

DOUBLE  PRECISION  A, B,C,D, E, EPS . XH , P, Q, S,TOL1 ,T0L2,U ,Y , W 

DOUBLE  PRECISION  ?U,FV,FW,FX,X 

DOUBLE  PRECISION  DABS,DSQBT,DSIGN 

COMMON  /SEARCH/  EPS 

C  *  0.  5D0*  (3.  -  DSQRT  (5.0D0)  ) 

A  *  AX 
B  *  BX 


00002510 

00002520 

00002530 

00002540 

00002550 

00002560 

00002570 

00002580 

00002590 

00002600 

00002610 

00002620 

00002630 

00002640 

00002650 

00002660 

00002670 

00002680 

00002690 

00002700 

00002710 

00002720 

00002730 

00002740 

00002750 

00002760 

00002770 

00002780 

00002790 

00002800 

00002810 

00002820 

00002830 

00002840 

00002850 

00002860 

00002870 

00002880 

00002890 

00002900 

00002910 

00002920 

00002930 

00002940 

00002950 

00002960 

00002970 

00002980 

00002990 

00003000 


7  =  A  ♦  C*  (B  -  A) 

W  =  7 
X  -  7 
2  =  O.ODO 
PX  *  F  (X) 

FV  =  FX 
FW  a  ?X 

20  XM  *  0.5D0*  (A  ♦  B) 

T0L1  *  EPS*DABS(X)  +  TOL/3.0D0 
T0L2  =  2. ODO*TOL1 

IF  (DABS  (X  -  XM)  .LB.  (T0L2  -  0.5D0*(B  -A)))  GO  TO  90 
IF  (DABS  {£)  .LE.  TOL1)  GO  TO  40 
B  ■  .  (X  -  W)  *(PX  -  F7) 

■  Q  =  (T  -  V)  *(F X  -  FW) 

P  =  (X  -  7)  *Q  -  (X  -  W)  *R 

Q  »  2.0D00*  (Q  -  R) 

IP  (Q  .GT.  O.ODO)  P  *  -P 
Q  =  DABS  (Q) 

R  =  E 
E  =  D 

30  IF  (DABS  (P)-  .GE.  DABS  (0.5D0*Q*R) )  GO  TO  40 
IF  (P  .LE.  Q*  (A  -X))  GO  TO  40 

IF  (P  .GE.  Q*(B  -  X))  GO  TO  40 

D  =  P/Q 
0  =  X  ♦  D 

IF  (  (0  -  A)  .IT.  TOL2)  D  *  DSIGN(TOL1,  XM  -  X) 

IF  (  (B  —  0)  .LT.  TOL2)  D  =  DSIGN(T0L1,  XM  -  X) 

GO  TO  50 

40  IF  (X  .GE.  XM)  E  =  A  -  X 

IF  (X  .LT.  XM)  E  =  B  -  X 

D  =  C*E 

50  IF  (DABS  (D)  .GE.  TOL1)  U  =  X  ♦  D 

IF  (DABS  (D)  .LT.  TOL1)  0  =  X  ♦  DSIGN(TOL1,  D) 

FU  *  P(TJ) 

IF  (FO  .GT.  FX)  GO  TO  €0 
IF  (0  .GE.  X)  A  =  X 

IF  (0  .LT.  X)  B  ■  X 

7  =  W 
F7  *  FW 
W  =  X 
FW  =  FX 
X  »  0 
FX  *  FO 
GO  TO  20 

60  IF  (0  .LT.  X)  A  *  0 

IF  (0  .GE.  X)  B  =  0 

IF  (FO  .LE.  FW)  GO  TO  70 
IF  (W  .EQ.  X)  GO  TO  70 
IF  (FO  .LE.  F7)  GO  TO  80 
IF  (7  .SQ.  X)  GO  TO  80 


00003010 

00003020 

00003030 

00003040 

00003050 

00003060 

00003070 

00003080 

00003090 

00003100 

00003110 

00003120 

00003130 

00003140 

00003150 

00003160 

00003170 

00003180 

00003190 

00003200 

00003210 

00003220 

00003230 

00003240 

00003250 

00003260 

00003270 

00003280 

00003290 

00003300 

00003310 

00003320 

00003330 

00003340 

00003350 

00003360 

00003370 

00003380 

00003390 

00003400 

00003410 

00003420 

00003430 

00003440 

00003450 

00003460 

00003470 

00003480 

00003490 

00003500 


IP  (7  .EQ.  H)  GO  TO  SO 

00003510 

GO  TO  20 

00003520 

70  7  =  W 

00003530 

-a 

< 

ii 

« 

00003540 

W  =  0 

00003550 

PH  *  FU 

00003560 

GO  TO  20 

00003570 

8  0  7  *0 

00003580 

?7  =  PO 

00003590 

GO  TO  20 

00003600 

90  PSIN  =  X 

00003610 

RETURN 

00003620 

END 

00003630 

;data- 

00003640 

5 

00003650 

3.85D0  -2.35D0  4.52D0 

2. 10600 

1.694D0 

00003660 

1.0D-68  1.0D0  1.0D1 

1. 0D-2 

1.0D3 

00003670 

3*********  ******************************************************** ******0000 00 10 

:  *00000020 

3  PROGRAM  NAME  :  SOLVALL  *00000030 

:  WRITTEN  BY  :  J.F.  BRANDEAU  *00000040 

3  COMPILER (S)  :  WATFIV  (DOOBLE  PRECISION)  *00000050 

:  *00000060 

:  PURPOSE  :  TO  CONVERT  EQUATION  STRAIN  =  F (RATE,  STRESS)  TO  EQUATION  *00000070 
3  STRESS  *  F (RATE,  STRAIN)  USING  SUBROUTINE  ZEROIN  TO  SOLVE  FOR  ROOT  OP*00000080 
E  THE  EQUATION.  CAN  HANDLE  UP  TO  SIX  SETS  OF  COEFFICIENTS  X  TO  PRODUCE* 00000090 


3  UP  TO  SIX  DATA  PAIRS  FOR  SAS  TO  GRAPH  ON  A  SINGLE  GRAPH. 


3  VARIABLES  : 


KGfi APS 


NUMBER  OF  SETS  OF  COEFFICIENTS  TO  BE  TAKEN  FROM  DATA. 


3  RATE  :  STRAIN  RATE  TO  BE  USED. 

3  KEPS  :  NUMBER  OF  POINTS  TO  BE  GENERATED  FCR  EACH  UNIQUE  VECTOR 
3  TOL  :  CONVERGENCE  CRITERION  FOR  ZEROIN. 

3  A,  B,  N,  D,  C  (I)  -  SET  OF  COEFFICIENTS  FCR  EACH  CURVE.  THE  I' 
3  ELEMENT  OF  EACH  CONSITITUTES  ONE  SET  OF  COEFFICIENTS  X. 


♦00000100 

*00000110 

*00000120 

*00000130 

*00000140 

*00000150 

*00000160 

*00000170 

*00000180 

*00000190 

*00000200 

*00000210 

*00000220 

*00000230 

*00000240 

*00000250 

*00000260 


FILE  #1 
FILE  #6 


3  SIG MAX  :  MAXIMUM  VALUE  OF  STRESS  (KSI)  EXEECTED  FOR  EACH  SET  OF  *00000250 

3  COEFFICIENTS  X.  THIS  IS  THE  HIGH  LIMIT  SENT  TO  ZEROIN.  THESE  HAY  *00000260 
3  BE  ADJUSTED  DOWNWARD  BY  THE  PROGRAM  IF  NEEDED  TO  PREVENT  UNDER/OVER  *00000270 

3  FLOWS.  IF  THE  VALUE  OF  SIGMIX(I)  IS  NOT  HIGH  ENOUGH  FOR  COEFFICIENT*00000280 

3  X(I),  THE  CURVE  WILL  BE  FLATTENED  AT  THE  HIGHER  STRESSES.  *00000290 

3  *00000300 

3  EPSHAX  :  MAXIMUM  VALUE  OF  STRAIN  FOE  WHICH  EACH  CURVE  IS  TO  BE  *00000310 

C  EVALUATED.  *00000320 

-  *00000330 

3  I/O  REQUIREMENTS  :  *00000340 

3  *00000350 

3  FILE  #1  :  ALL  INPUT  FROM  ABOVE,  FOLLOWING  SDATA  CARD.  *00000360 

3  FILE  #6  :  OUTPUT  OF  STR ESS- STRAIN  PAIRS  FCR  USE  BY  SAS  (LRSCL=130).  *00000370 

3  - t -  *00000380 

3  *00000390 

3  CCOMP  OPTIONS  (FORM  CSOPTIONS  CCCHP=?????)  :  *00000400 

3  *00000410 

3  4  :  OUTPUT  OF  RETURNED  STRESS  VALUES.  *00000420 

3  *00000430 

3*************************  *****************  ***********  ***************  **  *00000440 
3SOPTIONS  CCOHP*0  00000450 

IMPLICIT  REAL  *  8  (A-H,  N,  O-Z)  00000460 

EXTERNAL  FUNCT  00000470 

DIMENSION  A  (6)  ,  B  (6)  ,  C(6),  D(6),  N(6),  SIGMAX{6),  EPSMAX(6)  00000480 

DIMENSION  EPS  (50,6) ,  STRESS(50,6)  00000490 

COMMON  A,  B,  C,  D,  N,  TLOG ,  FATE,  KOUNT,  J,  EPS,  RATEB,  BATED  00000500 


ALL  INPUT  FROM  ABOVE,  FOLLOWING  SDATA  CARD. 

OUTPUT  OF  STRESS-STRAIN  PAIRS  FCR  USE  BY  SAS  (LRSCL=130). 


3  CCOMP  OPTIONS  (FORM  CSOPTIONS  CCCHP=?????)  : 
3  4  :  OUTPUT  OF  RETURNED  STRESS  VALUES. 


COSNON  /SOB2/  GEPS 

CALCULATE  SACHINE  EPSILON 

GEPS  =  1.0D0 
4  GEPS  =  GEPS/2.0D0 
TOL 1  =  1.000  ♦  GEPS 
IF  (TOL 1  .GT.  1.0D0)  GO  TO  4 

READ  (1,*)  KGRAFS,  RATE#  STEPS,  TOL 

TEPS  *  0.0D0 

DO  10  I  =  1,  KGRAFS 

READ  (1,*)  A  (I)  ,  B  (I)  ,  N  (I )  ,  D(I),  C(I),  SIGSAX(I),  SPSMAX(I) 

CHECK  FOR  POSSIBLE  OVERFLOW  FOR  HIGH  N'S 
ALTER  SIGNAX  DOWNWARD  IF  NECESSARY 

RATEB  a  DLOG 10 (RATE  **  B(I)) 

6  AX  a  N  (I)  *  DLOGIO  (SIGHAX(I))  +  RATEB 

IF  (AX.  GT.  7S.ODO)  THEN  DO 
SIGNAX(I)  =  SIGHAX(I)  -  0.5DO 
GO  TO  6 
ENDIF 

10  CONTINUE 

BEGIN  HAIN  LOOP 

DO  30  KOUNT  =  1,  KGRAFS 

STRESS  ( 1  ,KOUNT)  =  EPS(1,KOUNT)  =  0.0D0 

BX  =  SIGNAX  (KOUNT) 

DEPS  *  EPSH AX (KOUNT)  /  DFLOAT (KEPS-1) 

RATED  =  RATE  **  D (KOUNT) 

RATEB  =  RATE  **  B (KOUNT) 

WATCHING  FOR  VALUES  OF  N  THAT  WILL  CAUSE  OVERFLOW  OR  UNDERFLOW 

IF  (A(KOUNT)  .BQ.  0.0D0)  THEN  DO 
TLOG  *  -80.0D0 
ELSE  DO 

TLOG  a.  DLOGIO  (A  (KOUNT)) 

ENDIF 

INNER  LOOP  FOR  EACH  SOLUTION 
DO  20  J  *  2,  KBPS 

EPS(J, KOUNT)  =  DFLOAT  ( J-  1)  *  DEPS 
AX  a  STRESS  (J-1,  KOUNT) 

IF  (J  .EQ.  2)  AX  »  (  c (KOUNT)  *  EPS  (J, KOUNT)  *  RATED)  /  1.3D0 
STRESS  (J, KOUNT)  =  ZEROIN  (AX,  BX ,  FUNCT,  TOL) 

4  PRINT, 'FROH  ZEROIN,  STRESS’, J,»  =’, STRESS (J , KOUNT) 
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00000520 
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00000540 
00000550 
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00000580 
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000006 1 0 
00000620 
00000630 
00000640 
00000650 
00000660 
00000670 
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00000690 
00000700 
00000710 
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20  CONTINUE 
30  CONTINUE 

DO  40  I  =  1,  KEPS 

40  WHITE  (*,100)  (EPS  (I,J)  ,  STRFSS(I,J),  J  *  1,  KGB AFS) 

50  STOP 

100  FORMAT  ( 1 0  { 1PD1  3.  5)  ) 

END 

DOUBLE  PRECISION  FUNCTION  FUNCT (STRESS) 

IMPLICIT  REAL  *  8  (A-H,  N,  O-Z) 

DIMENSION  A  (6)  ,  B(6),  C  (6)  ,  D(6),  N(6),  EPS(50,6) 

COMMON  A,  B,  C,  D,  N,  TLOG ,  BATE,  KOUNT,  J,  EPS,  RATSB,  BATED 
IF  (STRESS.  EQ.  0.0D0)  THEN  DO 
FUNCT  »  -EPS (J, FOUNT) 

•  RETURN 
ENDIF 

TEMPI  =  STRESS  /  (C  (KOUNT)  *  RATED) 

IF  (STRESS  .GT.  1.0D-1)  GO  TO  10 

IF  ((N  (KOUNT)  *  DLOG10  (STRESS) )  .  LT.  -60.0D0)  THEN  DO 
HOLD  *  -25.0D0 
GO  TO  15 
ENDIF 

10  TEMP2  *  STRESS  **  N (FOUNT)  *  HATEB 
HOLD  *  DLOGIO  (TEHP2) 

15  IF  ((HOLD  ♦  TLOG)  .LT.  -17.0D0)  THEN  DO 
PUNCT  =  TEMPI  -  EPS  (J, FOUNT) 

ELSE  DO 

FUNCT  =  (TEMPI  ♦  A  (KOUNT)  *  TBHP2)  -  EPS (J, KOUNT) 

ENDIF 

>5  RETURN  ;  END 
.'SOPTIONS  NOLIST 

DOUBLE  PRECISION  FUNCTION  ZEBOIN  (AX,BX,F,TOL) 

DOUBLE  PRECISION  AX,BX,F,TCL 

DOUBLE  PRECISION  A ,B,C,D,E,EPS,FA, FB , FC,TOL1 , XM,P,Q,R,S 
DOUBLE  PRECISION  DABS,DSIGH 
COMMON  /SUB 2/  EPS 
A  a  AX 
B  a  BX 
FA  a  F(A) 

FB  a  F(B) 

20  C  =  A 
FC  a  FA 
D  a  B  -  A 
E  a  D 

30  IF  (DABS (FC)  .GE.  DABS  (FB) )  GO  TO  40 


A  * 

B 

B  = 

C 

C  * 

A 

FA 

a  FB 

FB 

a  FC 

FC 

*  FA 

00001010 
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00001120 

00001130 

00001140 
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00001170 

00001180 

00001190 

00001200 

00001210 

00001220 

00001230 

00001240 

00001250 

00001260 

00001270 

00001280 

00001290 

00001300 

00001310 

00001320 

00001330 

00001340 

00001350 

00001360 

00001370 

00001380 

00001390 
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00001410 
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00001470 

00001480 

00001490 

00001500 


o>  ut  V* 


40  TOL 1  =  2. 0D0*E?S*D A3S (B)  ♦  0.5D0*TOL 
XM  =  .  *5*  (C  -  3) 

IF  (DABS  (IS)  .LE.  TOL1)  GO  TO  90 
IF  (FB  .EQ.  0.0D0)  GO  TO  90 
IF  (DABS  (E)  .LT.  TOL1)  GO  TO  70 
IF  (DABS  (FA)  .LE.  DABS  (FB)  )  GO  TO  70 
IF  (A  .HE.  C)  GO  TO  50 
S  =  F3/FA 
P  =  2.0D0*XH*S 
Q  *  1.0D0  -  S 
GO  TO  60 
50  Q  =  FA/FC 
g  =  FB/FC 
•S  *  FB/FA 

P  =  S*  (2. 0D0*XB*Q*  (Q  -  R)  -  (B  -  A)*(H  -  1.0D0)  ) 

Q  *  (Q  -  1 . 0D0)  *  (B  -  1 . 0D0 )  *  (S  -  1.0D0) 

60  IF  (P  .GT.  O.ODO)  Q  =  -Q 
P  =  DABS  (P) 

IF  (  (2. 0D0*P)  .GE.  (3.0D0*X8*Q  -  DABS  (T0L1*Q)  )  )  GO  TO  70 
IF  (P  .GE.  DABS (0. 5D0*E*Q) )  GO  TO  70 
E  =  D 
D  =  P/Q 

GO  TO  80  . 

70  D  »  xa 
E  =  D 
80  A  =  B 
FA  =  FB 

IF  (DABS  (D)  .GT.  TOL1)  B  *  B  *  D 

IF  (DABS  (D)  .LE.  TOL1)  B  =  B  *  DSIGN(T0L1,  XB) 

FB  =  F  (B) 

IF  (  (FB*  (FC/DABS  (FC)  ) )  .GT.  O.ODO)  GO  TO  20 

GO  TO  30  ' 

90  ZEBOIN  =  B 
BSTUHN 
EHD 

DATA 

1.0D0  25  1.0D-7 

. 1 1 8 3D- 1 3  -3.74D-1  6.5342D0  6.707D-2  3.5514D3  50.0D0  .016D0 

36.777D-13  -4.1267D-1  7.6617D0  5.67034D-2  2.20498D3  25.0D0  .0060D0 

9.052D-13  -1.426D-1  7.632D0  6.287D-2  2.284D3  20.0  .0070D0 

3.709D-68  -2.3357D0  45.2472D0  1.7977D-2  1.694D3  45.0D0  .036D0 

).0D0  1.0D0  1.0D0  1.34946D-2  2.6354D2  15.0D0  .030D0 
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******************************* ********  ********** ************ ******* ***ooo 000 10 

*00000020 

PROGRAM  NAME  :  BAKSOLV  *00000030 

WRITTEN  BY  :  J.  ?.  BRAND E AO  *00000040 

COMPILES (S)  :  WATFIV  (DOUBLE  PRECISION)  *00000050 

PURPOSE  :  CONVERT  EQUATION  STRAIN  =  F (SATE,  STRESS)  TO  EQUATION  *00000070 

STRESS  =  F (RATE,  STRAIN)  BY  USING  SUBROUTINE  ZERCIN  TO  SOLVE  FOR  ROOT*00000080 


OF  EQUATION.  PROGRAM  CHECKS  RATE  TO  DETERMINE  EACH  CHANGE  AND 
PREVENT  UNDER  /  OVER  FLOW.  6  UNIQUE  RATES  ARE  ALLOWED  IN  THE  INPUT 
LIST  FROM  FILE  *4. 

VARIABLES  : 

SIG MAI  :  GREATER  THAN  MAXIMUM  VALUE  OF  STRESS  EXPECTED  FOR  EACH 
OBSERVED  PATF.  THIS  IS  THE  UPPER  LIMIT  FOR  ROOT  SEARCH,  AND  IS 
ADJUSTED  TO  PREVENT  OVER  /UNDER  PLOW.  IF  THIS  IS  TOO  LOW  THE  CURVE 
FOR  THAT  STRAIN  RATE  WILL  BE  FLATTENED  AT  THE  TOP. 

TOL  :  CONVERGENCE  CRITERION  FOR  ZEROIN. 

X  5  COEFF  :  MANTISSA  AND  EXPONENT  OF  COEFFICIENT  VECTOP.  PROGRAM 
COMBINES  BOTH  INTO  X. 


REQUIREMENTS 


FILE 


FILE 


OBSERVED  VALUES  OF  DATA  AS  USED  FOR  OTHER  PROGRAMS.  *00000270 
TITLE  MUST  BE  ON  FIRST  RECORD,  FOLLOWED  BY  ONE  OBSERVATION  *00000280 
PER  RECORD;  STRAIN  RATE,  STRAIN  AND  STRESS  (IN  ORDER).  *00000290 
X  S  COEFF.  X  IS  THE  MANTISSA  AND  COEFF  THE  EXPONENT  OF  *00000300 
THE  VARIABLE  COEFFICIENTS.  PRODUCT  X  *  COEFF  SHOULD  EQUAL  *00000310 
THE  COEFFICIENT.  *00000320 
OUTPUT  OF  POINTS  FOR  USE  BY  SAS.  *00000330 


*00000090 

*00000100 

*00000110 

*00000120 

♦00000130 

*-00000140 

*00000150 

*00000160 

*00000170 

*00000180 

*00000190 

*00000200 

*00000210 

*00000220 

*00000230 

*00000240 

*00000250 

*00000260 

*00000270 

*00000280 

*00000290 

*00000300 


FILE  #6  :  OUTPUT  OF  POINTS  FOR  USE  BY  SAS.  *00000330 

*00000340 

********************************* ********* **************************** *00 0003 50 
IMPLICIT  REAL  *  8  (A-H,  O-Z)  00000360 

EXTERNAL  FUNCT  00000370 

DIMENSION  X  (5)  ,  SIGMAX(6),  COEFF(5)  00000380 

COMMON  X,  HATE,  EPS,  RATEB ,  RATED  00000390 


INTEGER  TITLE (20) 

DATA  SIGMAX  /  50.0,  50 
TOL  =  1.0D-7 
READ  (5,*)  X 
READ  (5,*)  COEFF 
READ  (4,200)  TITLE 
WRITE  (3,300)  TITLE 
DO  5  J  *  1,  5 
X  (J)  =  X  (J)  *  COEFF  (J) 
WRITE  (3,100)  X(J) 
FATE1  =  0.0 


50.0, 


50.0,  50.0/ 


00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 


KP  *  0 

10  READ  (4,*,END=50)  RATE,  EPS,  STRESS 

:  BAS  STRAIN  RATE  CHANGED  IN  INPUT  LIST? 

IE  (RATS  .NE.  RATE1)  THEN  DO 
RATE1  *  RATE 
KP  =  KP  ♦  1 
KOUNT  =  0 

Z  CORRECT  SIGMAX  TO  PREVENT  OVER  /  UNDER  FLOW. 

RAT5B  =  DLOGIO  (RATE  **  X  f  2)  ) 

15- AX  *  X(3)  *  DLOGIO  (SIGMAX  (KP)}  ♦  RATEE 
IF  (AX  .GT.  75.0D0)  THEN  DO 
SIGMAX  (KP)  =  SIGMAX (KP)  -  0.5D0 
GO  TO  15 
ENDIP 

RATEB  =  BATE  **  X(2) 

RATED  =  RATE  **  X(4) 

ENDIF 

KOUNT  =  KOUNT  ♦  1 
IF  (EPS  .EO.  0.0D0)  THEN  DO 
SIG 1  *  0.0D0 
ELSE  DO 
AX  *  SIGl 
BX  »  SIGMAX  (KP) 

IF  (KOUNT  .EQ.  2)  AX  =  STRESS  /  1.3 
SIGl  =  ZEROIN  (  AX,  BX,  FUNCI,  TOL) 

ENDIF 

WHITE  (6,400)  HATE,  STRESS,  EPS,  SIGl 
GO  TO  10 
50  STOP 

100  FORMAT  (1H  ,1PD21.13) 

200  FORMAT  (20A4) 

300  FORMAT  (1H  , 'TITLE  :  ',20A4) 

400  FORMAT  (1H  ,5(1PD13.5)) 

END 

DOUBLE  PRECISION  FUNCTION  FUNCT ( STBES S) 

IMPLICIT  REAL  *  3  (A-H,  O-Z) 

DIMENSION  X  (5) 

COMMON  X,  RATE,  EPS,  RATEB,  RATED 
TEMPI  =  STRESS  /  X(5)  /  RATED 

FUNCT  =  TEMPI  ♦  X(1)  *  RATEB  *  STRESS  **  X  (3)  -  EPS 

RETURN 

END 

:$OPTIONS  NOLIST 

DOUBLE  PRECISION  FUNCTION  ZEROIN  (AX, B  X,P, TOL) 

DOUBLE  PRECISION  AX,BX,F,TOL 
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00000830 
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00000860 
00000870 
00000880 
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00000900 
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00000940 
00000950 
00000960 
00000970 
00000980 
00000990 
00001000 


DOOBLE  PRECISION  A , B, C, D, I , EPS , PA , FB , PC , TOL1 , XM , P , Q , P , S 
DOUBLE  PRECISION  EABS.DSIGN 
EPS  =  1.0D0 
10  EPS  =  EPS/2. 0D0 
TOL1  *  1.0D0  ♦  EPS 
IF  (TOL 1  .GT.  1.0D0)  GO  TO  10 
A  =  AI 
B  3  BE 
FA  3  F  (A) 

FB  -  ?(B) 

20  C  =  A 
FC  =  FA 
D  3  B  -  A 
•  E  =  D 

30  IF  (DABS (FC)  .GE.  DABS  (FB))  GO  TO  4 0 
A  =  B 
B  =  C 
C  3  A 
FA  =  FB 

FB  3  FC 

FC  =  FA 

40  TOL1  3  2. ODO*EPS*D ABS  (B)  ♦  0.5D0*TOL 

XW  *  .5*  (C  -  B) 

IF  (DABS  (XH)  .LE.  TOL1)  GO  TO  90 
IF  (FB  .EQ.  0.0D0)  GO  TO  90 
IF  (DABS  (E)  -IT.  T011)  GO  TO  70 
IF  (DABS (FA)  ;LE.  DABS  (FB) )  GO  TO  70 
IF  (A  .NE.  C)  GO  TO  50 
S  =  FB/FA 
P  =  2.0D0*XN*S 
Q  =  1.0D0  -  S 
GO  TO  60 
50  Q  =  FA/FC 
R  =  FB/FC 
S  =  FB/FA 

P  3  S*(2.0D0*XH*Q*(Q  -  E)  -  (B  -  A )  *  ( B  -  1.000)) 

Q  3  (Q  -  1 . 0D0)  *  (R  -  1 . 0D0 )  *  (S  -  1.0D0) 

60  IF  (P  .GT.  O.ODO)  Q  3  -Q 
P  3  DABS  (P) 

IF  (  (2.  0D0*P)  .GE.  ( 3. 0D0*XJ!*Q  -  DABS  (TOL1*Q)  )  )  GO  TO  70 
IF  (P  .GE.  DABS (0.5D0*E*Q) )  GO  TO  70 
E  3  D 
D  3  P/Q 
GO  TO  80 
70  D  *  XH 
B  3  D 
80  A  3  B 
FA  3  FB 

IF  (DABS  (D)  .GT.  TOL1)  B  3  B  ♦  D 

IF  (DABS  (D)  .LE.  TOL1)  B  3  B  ♦  DSIGN(TOL1,  XS) 
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PB  =  F (B) 

IP  (  (FB*  (FC/DABS  (FC)  )  )  .GT.  O.ODO) 
GO  TO  30 
90  ZEROIN  =  B 
RETtJRN 
END 


00001510 

GO  TO  20  00001520 

00001530 

00001540 

00001550 

00001560 


************  ************************************ ***** ********** ******* *000000 1 0 

*00000020 

PROGRAM  NAME  :  MAP2  *00000030 

WRITTEN  BY  :  J.F.  B»ANDEAU  *00000040 

COMPILER  (S)  :  WATFI7  (DOOBLE  PRECISION)  *00000050 


*00000030 
*00000040 
*00000050 
*00000060 

PORPOSS  :  USE  A  SET  OF  OBSERVED  POINTS  AND  COEFFICIENT  VECTOR  X  TO  *00000070 
HATCH  A  PREDICTED  VALUE  OF  STRAIN  =  F  ( STRESS ,  STRAIN  RATE)  TO  EACH  *00000080 
OBSERVED  POINT.  READS  FROM  EXTERNAL  FILE  FOR  COEFFICIENTS  AND  OBS-  *00000090 
SHVED  VALUES,  AND  CAN  READ  COEFFICIENTS  FROM  TRAILING  LIST  (THESE  *00000100 

OVERRIDE  EXTERNAL  VALUES).  SENDS  OBSERVED  POINTS  AND  PREDICTED  *00000110 

POINTS  TO  A  FILE  TO  BE  USED  BY  SAS.  CAN  ALSO  DO  SENSITIVITY  ANALYSIS*00000120 
WITHOUT  CHANGING  EXTERNAL  VALUES  OF  COEFFICIENTS.  *00000130 

VARIABLES  :  *00000150 

*00000160 

X  6  COEFF  :  MANTISSAS  AND  EXPONENTS  OF  COEFFICIENT  VECTOR.  THESE  ARE*00000170 
COMBINED  INTO  X  IN  THE  PROGRAM.  *00000190 

*00000190 

DELTA  :  FRACTIONAL  CHANGE  IN  VARIABLE  K  FOR  SENSITIVITY  ANALYSIS.  *00000200 
SET  THIS  EQUAL  TO  ZERO  TO  GET  TRUE  COEFFICIENTS.  *00000210 

*00000220 

K  :  VARIABLE  THAT  WILL  BE  ALTEHED  BY  AMOUNT  (DELTA  *  X(K)).  MUST  BE  *00000230 
BETWEEN  1  AND  5  ALWAYS.  *00000240 

*00000250 

STRESS  :  OBSERVED  VALUES  OF  STRESS  (KSI) .  *00000260 

*00000270 

HATE  J  OBSERVED  VALUES  OF  STRAIN  BATE  (1/SEC).  *00000280 

*00000290 

EPS  :  OBSERVED  VALUES  OF  STRAIN  (IN/IN).  *00000300 

*00000310 

EPS1  :  PREDICTED  VALUE  OF  STRAIN  RATE  (IN/IN)  .  *00000320 

-  *00000330 

I/O  REQUIREMENTS  :  *00000340 

*00000350 

FILE  #1  :  (OPTIONAL)  X  6  COEFF  VECTORS  ON  TWO  RECORDS.  *00000360 

FILE  *4  :  OBSERVED  VALUES  OF  RATE,  EPS,  STRESS  IN  THIS  ORDER.  THE  *00000370 
FIRST  CARD  MUST  BE  A  TITLE.  THE  REMAINING  CARDS  CONTAIN  *00000380 
ONE  OBSERVATION  EACH.  *00000390 

FILE  #5  :  X  6  COEFF  VECTORS  ON  TWO  RECORDS.  WILL  ALWAYS  BE  READ.  *00000400 

-  *00000410 

-  OPTIONS  (FORM  CSOPTIONS  CCOHP= ??????)  :  *00000420 

*00000430 

3  :  READ  FROM  DATA  CARDS  AT  END  OF  PROGRAM  LIST,  FOLLOWING  $DATA  CARD*00000440 

*00000450 

******************************************************************** ***Q0000Q60 
SOPTIONS  CCOHP=0  00000470 

IMPLICIT  REAL  *  8  (A-H,  N,  O-Z)  00000480 

DIMENSION  X  (5)  ,  COEFF  (5)  00000490 

INTEGER  TITLE  (20)  00000500 


DELTA  :  FRACTIONAL  CHANGE  IN  VARIABLE  K  FOR  SENSITIVITY  ANALYSIS. 
SET  THIS  EQUAL  TO  ZERO  TO  GET  TRUE  COEFFICIENTS. 


K  :  VARIABLE  THAT  WILL  BE  ALTEHED  BY  AMOUNT  (DELTA 
BETWEEN  1  AND  5  ALWAYS. 

STRESS  :  OBSERVED  VALUES  OF  STRESS  (KSI). 

HATE  :  OBSERVED  VALUES  OF  STRAIN  BATE  (1/SEC). 

EPS  :  OBSERVED  VALUES  OF  STRAIN  (IN/IN). 

EPS1  :  PREDICTED  VALUE  OF  STRAIN  RATE  (IN/IN). 

I/O  REQUIREMENTS  : 


X(K)) 


FILE  *1 
FILE  *4 


FILE  #5 


(OPTIONAL)  X  6  COEFF  VECTORS  ON  TWO  RECORDS. 

OBSERVED  VALUES  OF  RATE,  EPS,  STRESS  IN  THIS  ORDER.  THE 
FIRST  CARD  MUST  BE  A  TITLE.  THE  REMAINING  CARDS  CONTAIN 
ONE  OBSERVATION  EACH. 

X  6  COEFF  VECTORS  ON  TWO  RECORDS.  WILL  ALWAYS  BE  READ. 


-  OPTIONS  (FORM  CSOPTIONS  CCOHP= ??????) 


READ  (5,*)  X 
READ  (5,*)  COEFF 
:3  READ  (1,*)  X 

:3  READ  (1,*)  COEFF 

DELTA  =  O.ODO 
K  -  1 

X{K)  m  X(K)  ♦  (1.0D0  +  DELTA) 

DO  5  J  s  1,  5 

X(J)  *  X(J)  *  COEFF  (J) 

5  WHITE  (3,75)  X ( J) 

HEAD  (4,200)  TITLE 
WRITE  (3,300)  TITLE 

10  READ  (4,*,END=50)  FATE,  EPS,  STRESS 

:  CALCULATE  PREDIICTED  STRAIN 

TEMPI  *  STRESS  /  X  (5)  /  RATE  **  X  ( 4) 
TE3P2  =  STRESS  **  X(3)  *  RATE  **  X(2) 
EPS1  =  TEMPI  +  X(1)  *  TE8P2 

WRITE  (6,1Q0)  RATE,  STRESS,  EPS,  EPS1 
GO  TO  10 
50  STOP 

75  FORMAT  (1H  ,  1PD20.12) 

100  FORMAT  (4D25.  13) 

200  FORMAT  f20A4) 

300  FORMAT  (1H  ,20A4) 

END 

5DATA 

5. 63  -0.  589  6.19  3.20  3.10 

1.0D-13  1.0D0  1.0D0  1.0D-2  1.0D3 


00000510 

00000520 

00000530 

00000540 

00000550 

00000560 

00000570 

00000580 

00000590 

00000600 

00000610 

00000620 

00000630 

00000640 

00000650 

00000660 

00000670 

00000680 

00000690 

00000700 

00000710 

OOU00720 

00000730 

00000740 

00000750 

00000760 

00000770 

00000780 

00000790 

00000800 

00000810 


/F0PCE1  JOB  DU.D08.  AQ0221  ,BP.ANDEAU,T=  (,  10)  ,M=  (2,0)  00000010 

/  EXEC  SATFIV  00000020 

/GO.  FT09F00  1  DD  DSN=DU. DOB. AQ0221 . BRAND EAU. DATA. ONE, DISP=SHP  00000030 

/GO. FT08F001  DD  DSN=DU. DOB. AQO 22 1. BRAND EAU. ATBDATA. ONE, DISP=SHR  00000040 

/GO.SYSIN  DD  *  00000050 

JOB  00000060 

SOPTIONS  CCOHP=2,NOEXT, NOCHECK, DEC  00000070 

***************  ***************  *******************  ************4>*****  ***Q00  00080 
SUBROUTINE  FORCE  IN  MAIN  PROGRAM  FORM  00000090 

THIS  PPOG5AM  ANALYZES  THE  ATB  OOTPOT  DATA  FOB  DATA  NEEDED  IN  THE  00000100 

OPERATION  OF  THE  BREAK  PROGRAM.  00000110 

1)  THE  TOP  OF  EACH  TIME  INCREMENT  DATA  SET  IS  FOUND  00000120 

2)  ROTATIONS  AND  ANGULAR  VELOCITIES  FOR  TIME  STEP  ARE  FOUND  00000130 

3)  DISPLACEMENTS  AND  LINEAR  VELOCITIES  ARE  FOUND  00000140 

4)  JOINT  FORCES  AND  TORQUES  ARE  FOUND  00000150 

5)  DATA  IS  WHITTEN  TO  PRINTER  AND/OR  DISK  00000160 

A.  OUTPUT  TO  DISK  IS  IN  UNFORMATTED  FORM  00000170 

B.  FOP.  EACH  TIME  INCREMENT,  THE  TIME  (MS)  AND  DATA  FOR  EACH  00000180 

LIMB  IS  OUTPUT.  EACH  IIMB  00000190 

DATA  SET  FOR  THAT  TIME  INCREMENT  IS  ON  A  RECORD,  PRECEDED  00000200 
BY  THE  IDENTIFYING  NUMBER  FOR  THAT  LIMB  (1  -  8).  00000210 

**********************************************************************00000220 
KCON  IS  THE  NUMBER  OF  "OTHER  CONSTRAINT  FORCES"  00000230 

-1  =  NONE  00000240 

>0  =  NUMBER  OF  ROWS  OF  DATA  TO  BE  FOUND  FOR  EACH  TIME  STEP  00000250 

00000260 

POSIT  =  POINT  OF  ATTACHMENT  RELATIVE  TO  C.G.  OF  SEGMENT  IN  SEGMENT  00000270 

LOCAL  Z- AXIAL  COORDINATES.  00000280 

CONSTRN  (I, J)  =  FORCES  IN  INERTIAL  COORDINATES  (J  =  1,  3)  FOR  SEGMENT  100000290 
(1*1,  KCON)  00000300 

00000310 

CHARACTER  *4  IDUM,IFLAG  00000320 

REAL  D  (8 ,30)  ,  VEH(18),  POSIT(24),  CNSTBN  (8,3)  00000330 

REAL  *  8  DD  (3)  ,  DA  (3)  00000340 

INTEGER  KT3  (8)  00000350 

DATA  KT3  /  8  *  1  /,  CNSTBN  /24  *  0.0/  00000360 

NT  =  31  ;  DLT  =  0.01  00000370 

IW  »  3  ;  IR  =  1  ;  IDISK1  =  9  ;  IDISK2  =  8  00000380 

WHITE  (IW,1*00)  00000390 

JP  a  1  00000400 

READ  (1,*)  KCON  00000410 

KCON 4  *  KCON  *  4  00000420 

IF  (KCON4  .EQ.  0)  KCON4  *  1  00000430 

POSIT (1)  »  0  00000440 

DO  2  I  *  1,  KCON 4,  4  00000450 

READ  (1,*)  POSIT  (I)  ,  POSIT  (1+1)  ,  POSIT(I+2),  POSIT  (1+3)  00000460 

2  KT3  (POSIT  (I))  »  3  00000470 

.2  WRITE  (IDISK2)  NT,  DLT,  KC0N4,  (POSIT  (I),  1=1,  KCOH4)  00000480 

DO  4  J  =  1,  200  ”  00000490 

READ  (IDISK1 , 400)  IDUH  00000500 


3 


2 


IF  (IDUH  . EQ.  *  5  H')  GO  TC  6 

00000510 

4  CONTINUE 

00000520 

6  DO  8  J  =  1,  10 

00000530 

READ  (IDI SKI , 420)  I,  W,  I,  Y,  Z,  A,  B,  C 

00000540 

IF  (J  .EQ.  3  .OR.  J  .EQ  .6)  GO  TO  8 

00000550 

WHITE  {IW,  120)  I,  W,  X,  Y,  Z,  Af  B,  C 

00000560 

00000570 

OUTPUT  TO  DISK  IN  X-AXIAL  COORDINATES  AND  CONSECUTIVE  SEGHENT 

NUHBERS  00000580 
00000590 

WRITE  (IDISK2)  JP,  W,  Z,  X,  Y ,  C,  A,  B 

00000600 

JP  =  JP  ♦  1 

00000610 

8  CONTINUE 

00000620 

WRITE  (IW,  130)  (POSIT  (KK)  ,  KK  =  1,  KCON4) 

00000630 

•  IP  =  0 

00000640 

DO  80  I  =  1,  1000 

00000650 

ON  ERROR  GOTO  75 

00000660 

IF  (IP  .GE.  NT)  THEN  DO 

00000670 

PRINT, IP,'  TIBS  STEPS  FOUND* 

00000680 

STOP  ;  ENDIF 

00000690 

00000700 

FIND  TOP  OF  DATA  SET  FOR  EACH  TIHE  INCREHENT 

00000710 

00000720 

13  READ  (IDISK1  ,900,  END  *  90)  IFLAG 

00000730 

IF  (IFLAG  .  NE.  '  HAIN')  GO  TO  13 

00000740 

BACKSPACE  IDISK1 

00000750 

READ  (IDISK1,  905)  TIHE 

00000760 

TIHE  =  TIHE  /  1000. 

00000770 

00000780 

FIND  ROTATIONS  ,  ANGULAR  VEL.  AND  ANGULAR  ACC.  FOR  THIS  TIHE 

STEP 

00000790 

00000800 

IP  *  IP  ♦  1 

00000810 

DO  20  II  =  1,  15 

00000820 

HEAD  (IDISK1,  1000)  IDUH 

00000830 

IF  (IDUH  .NE.  'H  •)  GO  TO  20 

00000840 

JP  *  1 

00000850 

DO  15  J  *  1,  10 

00000860 

READ  (IDxSKI  ,1100)  (D(JP,KK),KK  =  7,9),  (DD  (KK)  ,  KK  =  1, 

3)  , 

00000870 

*  (DA  (KK)  ,  KK  =  1,3) 

00000880 

IF  (J  .EQ.  3  .OR.  J  .EQ.  6)  GO  TO  15 

00000890 

DO  1»  K#  «  1,  3 

00000900 

D (JP, K4 *18)  =  SNGL  (DA  (K4) ) 

00000910 

14  D(JP,K4+9)  =  SNGL  (DD  (K4) ) 

00000920 

JP  *  JP  ♦  1 

00000930 

15  CONTINUE 

00000940 

HEAD  (IDISK  1,800)  (VEH(KK),  KK  =  7 ,  12) 

00000950 

GO  TO  22 

00000960 

20  CONTINUE 

00000970 

- 

00000980 

FIND  LINEAR  POSITION,  VELOCITY  AND  ACCELERATION  FOR  THIS  TIHE 

STEP 

00000990 

00001000 

22  DO  UO  II  =  1,  15 

HEAD  (IDISK  1  ,1000)  IDUH 
IF  (IDUH  .ME.  * H  •)  GO  TO  40 

JP  =  1 

DO  25  J  =  1,  10 

HEAD  (IDISK  1,200  0)  (D(JP,KK),KK  =  1,  6),  (D(JP,KK),  KK  -  22, 
IF  (  J  -EQ.  3  .OR.  J  .  EQ.  6  )  GO  TO  25 
JP  =  JP  +  1 
25  CONTINUE 

HEAD  (IDISK  1, 1000)  IDUH 

HEAD  (IDISK  1, 2000)  (VEH(KK),  KK  =  1,  6) 

GO  TO  50 
40  CONTINUE 

FIND  U1  5  U2  ARRAYS 

50  DO  42  II  =  1,  25 

HEAD  (IDI SKI , 400)  IDUH 
42  IF  (IDUH  .EQ.  *  5  H')  GO  TO  44 
44  JP  =  1 

DO  46  JJ  =  -1,  10 

READ  (IDISK1 , 2000)  (D(JP,KK),  KK  =  25,  30) 

I?  (JJ  .EQ.  3  .OR.  JJ  .EQ.  6)  GO  TO  46 
JP  =  JP  ♦  1 
46  CONTINUE 

FIND  FORCES  AND  TORQUES  FOR  THIS  TIME  INCBEHENT 

t 

DO  70  IJ  =  1,  50 

READ  (IDI SKI  ,1000)  IDDK 
IF  (IDUH  .NS.  'HP  *)  GO  TO  70 
JP  =  1 

DO  60  JJ  »  1,  10 

READ  (IDISK1  ,4000)  (D(JP,KK),  KK  =  13,  18) 

IF  (JJ  .EQ.  3  .OR.  JJ  .EQ.  6)  GO  TO  60 
JP  =  JP  ♦  1 
60  CONTINUE 

READ  (IDISK1, 4000)  (VEH(KK)  ,  KK  *  13,  18) 

GO  TO  71 

70  CONTINUE 

71  DO  74  II  =  1,  20 

HEAD  (IDISK1 ,400)  IDUH 
IF  (IDUH  .HE.  •  NO.*)  GO  TO  74 
HEAD  (IDISK1 , 400)  IDUH 
DO  72  IJ  *  1,  KCON 

72  READ  (IDISK1,  500)  (CNSTBN  (POSIT  (4*IJ-3)  ,IN)  ,  IN  =  1 ,  3) 

GO  TO  75 

74  CONTINUE 

75  CONTINUE 


00001010 
00001020 
00001030 
00001040 
00001050 
24)  00001060 
00001070 
00001080 
00001090 
00001100 
00001110 
00001120 
00001130 
00001140 
00001150 
00001160 
00001170 
00001180 
00001190 
00001200 
00001210 
00001220 
00001230 
00001240 
00001250 
00001260 
00001270 
00001280 
00001290 
00001300 
00001310 
00001320 
00001330 
00001340 
00001350 
00001360 
00001370 
00001380 
00001390 
00001400 
00001410 
00001420 
00001430 
00001440 
00001450 
00001460 
00001470 
00001480 
00001490 
00001500 


:  OUTPUT  RESULTS 

:  C2  =  DISK  OUTPUT  (UNIT  =  IDISK2) 

3  C3  *  LINE  PRINTER  (UNIT  =  IN) 

12  WRITE  (IDISK2)  TIME,  7EH 

:3  PRINT, *  TIME  =',TIME,'  VEHICLE  =',  VEH 

32  PRINT, 'TIME  =',TIMB, '  SEC* 

DO  85  J=  1*8 
K3  =  ET3  ( J) 

22  WRITE  (IDISK2)  J,K3,  (D  (J,KK)  ,  KK  =  1 , 30)  ,  (CNSTRN  (J,  KK)  ,  KK=1 ,  K3) 

33  WRITE  (IW , 808)  J  ,K3,  (D  ( J  ,KK)  *  KK=1 , 30)  ,  (CNSTRN  ( J,  *K)  ,  KK=1 ,  K3) 
85  CONTINUE 

80  CONTINUE 

90'PRINT, 'END  OP  FILE  REACHED  ON  UNIT*,IDISK1 
PRINT, IP,*  TIME  STEPS  POUND' 

STOP 

3  INPUT  FORMATS 

130  FORMAT  ('  EXTRA  FORCE  LOCATIONS  :»,F5.0,  3F12.3) 

400  FORMAT  (IX, A4) 

500  FORMAT  (25X-,  3F15.5) 

420  FORMAT  (13,  131,  F7.3,  3X,  3F11.5,  5X ,  3F8.3) 

800  FORMAT  (/,  11X,  3F9.4,  41,  3  (F7.3,  71)) 

900  FORMAT  (7X,  A4) 

905  FORMAT  (57X,  F8.3) 

1000  FORMAT  ( 4X *  A4) 

1100  FORMAT  (11X,  3F9-4,  3X,  3D14.5,  31,  3D14.5) 

2000  FORMAT  (11X,  3F11.4,  3X,  3F12.5,  31,  3F14.5) 

4000  FORMAT  (15X,  3F11.4,  3X,  3F12.5) 

3  OUTPUT  FORMATS 

100  FORMAT  ( 1 H 1 ,  10X,  'FORCES  FROM  ATB  MODEL*/) 

33120  FORMAT  (*  I*',  12,  '  W=* ,  F7.3,  •  XTZ=',  3F10.5,  »  ABC=»,  3F7. 
33808  FORMAT  (1H  ,  12,  13,  2X,  (T1 1,6 {E '6.7,  3X) ) ) 

END 

SDATA 

2 

2  0.0  0.0  0.0 

4  0.0  0.0  0.0 

3JSTOP 

38  END 
/* 

V*PW*BOWE 

V 


00001510 
00001520 
00001530 
00001540 
00001550 
00001560 
00001570 
00001580 
00001590 
00001600 
00001610 
00001620 
00001630 
00001640 
00001650 
00001660 
00001670 
•00001680 
00001690 
00001700 
00001710 
00001720 
00001730 
00001740 
00001750 
00001760 
00001770 
00001780 
00001790 
00001800 
00001810 
00001820 
.3)  00001830 
00001840 
00001850 
00001860 
00001870 
00001880 
00001890 
00001900 
00001910 
00001920 
00001930 
00001940 


'/F0RCE2  JOB  DO.D08.AQ0221,BBAMDEAO,T=(,10)  ,P=45,M=  (2,0) 

V  EXEC  WATFIV 

'/GO.  FT09F00 1  DD  DSN=DU.  DOB.  AQO  221 .  BP. AND E  AO.  DATA. TWO  ,  DISP=SH5 
'/GO. FT08F001  DD  DSN=DU. D08. AQO 2 21 . BRAND E AD. ATBDATA. TWO, DISP=  SHE 
'/GO.SYSIN  DD  * 

5JOB 

:$OPTIONS  NOSXT,CCOMP=2, NOCHECK, DEC 
:$OPTIONS  NOLIST 

CHARACTER  *  10  IDUM,  DUMB 
CHARACTER  *  3  ISEG1,  ISEG2 

REAL  PANEL  (32,8,3,*),  SEG NET  (32,8,3,6),  XI  (3),  X2  (3) 

INTEGER  CHEKP  (32,8),  CHEKS  (32,8),  ID{3),  OCOONT,  NP(4) 

CHARACTER  ITEST*3 (3) /• RUL’ ,»  FLL*.  * LUL1, • LLL* , 1 RUA’ , ’RLA* , ’LUA* 

*  ,*LLAV 

DATA  CHEKP, CHEKS  /  512  *  0/ 

DATA  PANEL, SEGMET  /9216  *  0.0/ 

DATA  NP  /  58,  92,  92,  92  / 

********************************************-************************  ** 

EXPLANATION  OF  VARIABLES  -- 

PANEL  (I,J,K,L)  AND  SEGMET  (I,J,K,L)  CONTAIN  CONTACT  DATA 
I  =  TIME  STEP 
J  =  LIMB  OR  MEMBER  NUMBER 
K  =  CONTACT  NUMBER 
L  =  1  TO  6 

L  =  1  IS  THE  NUMBER  OF  THE  PANEL  OR  SEGMENT  CONTACTED 

—  SEGMENT  NUMBERS  ARE  IN  ATS  MODEL  CODE,  NOT  BREAK  CODE 
L  =  2  IS  THE  NORMAL  FORCE 
L  ■  3  IS  THE  FRICTION  FORCE 

L  =  #  TO  6  ARE  THE  X,Y,Z  COORDINATES  OF  THE  CONTACT  POINT 

SEGMENT  CONTACT  POINTS  ARE  IN  AT B  LOCAL  COORDINATES 
PANEL  CONTACT  POINTS  ARE  IN  ATB  INERTIAL  COORDINATES 

CHEKP  (I,J)  AND  CHEKS  (I,J)  ARE  CONTACT  COUNTERS  FOR  EACH  LIMB 
AND  TIME  STEP 
I  =  TIME  STEP 
J  =  LIMB  OR  MEMBER  NUMBER 

THE  VALUE  STORED  IN  LOCATION  CHEKP  OR  CHEKS  (I,J)  IS  THE  NUMBER 
OF  CONTACTS  FOR  THAT  LIMB  AND  TIME  STEP  THAT  HAD  NON-ZERO 
FORCES.  CONTACTS  WITH  ZERO  FORCES  ARE  NOT  STORED  IN  THE  ARRAY 
OR  WRITTEN  TO  THE  DISK  OR  PRINTER. 

NP  =  THE  NUMBER  0?  RECORDS  THAT  OCCUR  BETWEEN  PAGE  HEADERS 
ON  THE  ATB  OUTPUT  FILE.  THIS  MUST  BE  CHANGED  AS  THE 
TIME  STEP  THEY  USE  CHANGES.  THIS  IS  NOT  THE  SAME  BETWEEN 
EACH  HEADER  BLOCK. 

UCOUNT  *  THE  NUMBER  OF  CONTACT  FILES  TO  BE  LOOKED  FOR.  THE 
HUMBER  OF  FILES  FOUND  IS  COUNTED,  NOT  THE  NUMBER  OF 
INDIVIDUAL  CONTACTS.  UCOUNT  DOES  NOT  INCLUDE  FILES  HAVING 


00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

*00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250 

00000260 

00000270 

00000280 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 


ALL  ZERO  FORCES  OB  SEGMENTS  THAT  ARE  NOT  OF  INTEREST.  00000510 

00000520 

TMAX  =  MAXIMUM  TIMS  TO  BE  FOUND.  00000530 

00000540 

PORSAT  #  3416  MUST  ALSO  BE  CHANGED  TO  SNIP  THE  PROPER  NUMBER  OF  00000550 

RECORDS  IN  EACH  FILE  IF  THE  FILE  IS  CF  NO  INTEREST.  00000560 

00000570 

* ***** ********** ******* ************************************************00000580 


TMAX  =  MAXIMUM  TIMS  TO  BE  FOUND. 

PORMAT  #  3416  MUST  ALSO  BE  CHANGED  TO  SNIP  THE  PROPER  NUMBER  OF 
RECORDS  IN  EACH  FILE  IF  THE  FILE  IS  CF  NO  INTEREST. 


ON  ERROR  GOTO  255 

DLT  =  0.01  ;  NT  *  31;  THAX  *  300.0 

IW  *  3  IR  =  1  ;  IDISKI  =  9  ; 

WRITE  (IB,  100) 

IP  =  0  ;  IDT  =  INT  (10000.  *  DLT)  ; 

DO  200  1=1,  500 
READ  (IDISKI,  103,  END=255)  IDUM 
IF  (IDUM  .EQ.  *  SEGMENT  NO*  )  GO  TO  210 
IF  (IDUM  .NE.  ‘VEHICLE  PA»)  GO  TO  200 


00000590 

.0  00000600 
;  I DISK 2  =  8  00000610 

00000620 

;  ICOUNT  =  0  ;  UCOUNT  =  9  00000630 


READ  (IDISKI,  102)  IPAN1,  ISEG1,  IPAN2,  ISEG2 
IFLAG1  =  IFLAG2  =  0 

CHECK  FOR  SEGMENT  OF  INTEREST 


DO  200  1=1,  500  00000640 

READ  (IDISKI,  103,  END=255)  IDUM  00000650 

IF  (IDUM  .EQ.  ‘SEGMENT  NO*  )  GO  TO  210  00000660 

IF  (IDUM  .NE.  ‘VEHICLE  PA»)  GO  TO  200  00000670 

:  00000680 

I  -  PANEL  VS.  SEGMENT  CONTACT  -  00000690 

:  00000700 

READ  {IDISKi,  102)  IPAN1,  ISEG1,  IPAN2,  ISEG2  00000710 

IFLAG1  =  IFLAG2  =  0  00000720 

:  00000730 

:  CHECK  FOR  SEGMENT  OF  INTEREST  00000740 

I  00000750 

DO  110  J  «  1,  8  .  00000760 

IF  (ISEG1  .NE.  ITEST(J)  .AND.  IS2G2  .  NE.  ITEST  (J) )  GO  TO  110  00000770 
IF  (ISSG1  .NE.  ITEST  (J))  GOTO  105  00000780 

JHOLD1  =  J  ;  IFLAG1  =  1  00000790 

105  IF  (ISEG2  .NE.  ITEST(J))  GO  TO  109  00000800 

JHOLD2  =  J  ;  IFLAG2  =  2  00000810 

109  IFLAG3  =  I FLAG 1  ♦  IFLAG2  00000820 

IF  (IFLAG3  -  3)  110,  112,  110  00000830 

110  CONTINUE  *  00000840 

:  00000850 

IF  (IFLAG1  ♦  IFLAG2)  1T1,  111,  112  00000860 

:  00000870 

:  SKIP  ENTIRE  FILE  00000880 

:  00000890 

111  READ  (IDISKI, 3416)  DUMB  00000900 

GO  TO  200  00000910 

112  READ  (IDISKI, 113)  DOMB  00000920 

NOHPAG  =  NP(1)  .  00000930 

ICOUNT  =  ICOUNT  *  1  00000940 

DO  155  J  -  1,11  00000950 

IF  (J  .EQ.  1)  GO  TO  120  00000960 

READ  (IDISKI,  114)  DUMB  000009^0 

120  DO  150  JL=  1,  NUMPAG  00000980 

READ  (IDISKI,  115)  T,  FN1,  FF 1,  XI,  FN2,  FF2,  X2  00000990 

IT  =  INT  (10.  *  T)  00001000 


ITT  =  IT  /  IDT  *  1 

00001010 

IF  ((IT/IDT  *  IDT)  .NE.  IT)  GO  TO  150 

00001020 

IF  (IFLAG2)  500,  130,  129 

00001030 

29 

I?  (FN2  .EQ.  0.0  .AND.  FF 2  .EQ.  0.0)  GO  TO  131 

00001040 

CHEKP  {ITT,  JHOLD2)  =  TEMPCT  =  CHEKP  (ITT,  JHOLD2) 

*  1 

00001050 

PANEL  (ITT,  JHOLD2,  TEMPCT,  1)  =  FLOAT  (IPAN2) 

00001060 

PANEL  {ITT,  JHOLD2,  TEMPCT,  2)  =  7N2 

00001070 

PANEL  (ITT,  JHOLD2,  TEMPCT,  3)  =  FF2 

00001080 

DO  132  K4  =  1,  3 

00001090 

32 

PANEL  (ITT,  JHOLD2,  TEMPCT,  K4  +  3)  =  X2  (K4) 

00001100 

31 

I?  (IFLAG1)  500  ,  140  ,  130 

00001110 

30 

IF  (FN 1  .EQ.  0.0  .AND.  FF 1  .EQ.  0.0)  GO  TO  140 

00001120 

CHEKP  (ITT,  JHOLD1)  =  TEMPCT  =  CHEKP (ITT,  JHOLD1)  ♦ 

1 

00001130 

PANEL  (ITT,  JHOLD1,  TEMPCT,  1)  *  FLOAT  (IPAN1) 

00001140 

PANEL  (ITT,  JHOLD1 ,  TEMPCT,  2)  *  FN 1 

00001150 

PANEL  (ITT,  JHOLD1,  TEMPCT,  3)  =  FF 1 

00001160 

DO  133  K4  =  1,  3 

00001170 

133 

PANEL  (ITT,  JHOLD1 ,  TEMPCT,  K4+3)  =  XI  (K4) 

00001180 

140 

IF  (ITT  .GE.  NT  .OR.  T  „EQ.  TMAX)  GO  TO  200 

00001190 

150 

CONTINOE 

00001200 

NDMPAG  =  NP(J-H) 

00001210 

55 

CONTINOE 

00001220 

!00 

CONTINOE 

00001230 

• 

00001240 

:05 

READ  (IDISK 1 ,  103,  END  =  255)  IDUM 

00001250 

IF  (IDO H  .NE.  'SEGMENT  NO')  GO  TO  260 

00001260 

• 

00001270 

-  SEGMENT  VS.  SEGMENT  CONTACT  - 

00001280 

• 

00001290 

!10 

BACKSPACE  IDISK1 

00001300 

READ  (IDISK  1,  101)  IDOM,  IHOLD1,  ISEG1,  IHOLD2,  ISEG2 

00001310 

IFLAG1  =  IFLAG2  =  0 

00001320 

JHOLD1  =  JHOLD2  =  0 

00001330 

* 

00001340 

:  CHECK  FOP  SEGMENT (S)  OF  INTEREST 

00001350 

• 

00001360 

DO  220  J  =  1,  8 

00001370 

IF  (ISEG1  .NE.  ITEST(J)  .AND.  ISEG2  .  NE.  ITEST  (J) ) 

GO  TO 

220  00001380 

IF  (ISEG1  .NE.  ITEST  (J))  GO  TO  218 

00001390 

JHOLD1  =  J  ;  IFL1G1  =  1 

> 

00001400 

GO  TO  220 

00001410 

!18 

IF  (ISEG2  .NE.  ITEST  ( J)  )  GO  TO  219 

00001420 

JHOLD2  =  J  ;  IFLAG2  =  2 

00001430 

!19 

IFLAG3  *  IFLAG1  ♦  IFLAG2 

0000144Q 

IF  (IFLAG3  -  3)  220,  222,  220 

00001450 

:20 

CONTINOE 

00001460 

• 

00001470 

IF  (IFLAG2  ♦  IFLAG1 )  500,  221,  222 

1 

00001480 

• 

00001490 

• 

SKIP  ENTIRE  FILE 

00001500 

• 

00001510 

121 

READ  (IDISK1,3416)  DUMB 

00001520 

GO  TO  260 

00001530 

!22 

READ  (IDISK1 , 113)  DUMB 

00001540 

NUMPAG  =  NP  (1) 

00001550 

ICOUNT  =  ICOUNT  ♦  1 

00001560 

ITT  =  1 

00001570 

DO  252  J  =  1,11 

00001580 

IF  (J  .EQ.  1)  GO  TO  224 

00001590 

READ  (IDISK 1 , 223)  DUHE 

00001600 

!24 

DO  250  JL  =  1,  NUMPAG 

00001610 

READ  (IDISK  1 ,  225,  END=255)  T, 

FN1 ,  FF1,  XI,  12 

00001620 

IF  (FN 1  .EQ.  0.0  .AND.  FF1  .EQ. 

0.0)  GO  TO  240 

00001630 

IT  =  I NT  (  10.  *  T  ) 

00001640 

ITT  =  IT  /  IDT  1 

00001650 

IF  ((IT  /  IDT  *  IDT)  .NE.  IT)  GO 

TO  240 

00001660 

IF  (IFIAG2)  500,  230,  226 

00001670 

?  2  6 

CHEKS  (ITT,  JHOLD2)  =  TEHPCT  = 

CHEKS  (ITT,  JHOLD2)  ♦  1 

00001630 

SEGMET  (ITT,  JHOLD2,  TEMPCT,  1) 

*  FLOAT  (JHOLDl) 

00001690 

SEGMET  (ITT,  JHOLD2,  TEMPCT,  2) 

=  FN1 

00001700 

SEGMET  (ITT,  JHOLD2,  TEMPCT,  3) 

=  FF1 

00001710 

DO  227  K4  =  1,  3 

00001720 

’27 

SEGMET  (ITT,  JHOID2,  TEMPCT,  K4+3)  =  X2  (K4) 

00001730 

IF  (IF  LAG  1)  500,  240,  230 

00001740 

!30 

CHEKS  (ITT,  JHOLD1)  *  TEMPCT  = 

CHEKS  (ITT,  JHOLDl)  ♦  1 

00001750 

SEGMET  (ITT,  JHOLD1,  TEMPCT,  1) 

*  FLOAT  (JHOLD2) 

00001760 

SEGMET  (ITT,  JHOLD1 ,  TEMPCT,  2) 

*  FN1 

00001770 

SEGMET  (ITT,  JHOLDl,  TEMPCT,  3) 

=  FF1 

00001780 

DO  232  K 4  -  1 ,  3 

00001790 

?32 

SEGMET  (ITT,  JHOLDl,  TEMPCT,  K4  +  3)  =  XI  (K4) 

00001800 

240 

IF  (ITT.  GE.  NT  .OR.  T  .EQ.  T MAX)  GO  TO  260 

00001810 

150 

CONTINUE 

00001820 

NUMPAG  *  NP  (J*  1) 

00001830 

’52 

CONTINUE 

00001840 

00001850 

■60 

IF  (ICOUNT  -LT.  UCOUNT)  GO  TO  205 

00001860 

00001870 

:  OUTPUT  DATA 

000C1880 

• 

00001890 

?55 

DO  1000  J3  =  1,  NT 

00001900 

TIME  =  FLOAT  (J3-1)  *  DLT 

00001910 

12 

WRITE  (3,2200)  TIME 

00001920 

:3 

WRITE  (IDISK2)  TIME 

00001930 

DO  1000  J4  «  1,  8 

00001940 

ID  (1)  =  J4 

00001950 

IF  (CHEKP  (J3, J4)  .EQ.  0)  THEN  DO 

00001960 

JPAN  =  ID2  =  ID  (2)  = 

1 

00001970 

ELSE  DO 

00001980 

JPAN  =  6  ;  ID(2)  =  6 

*  CHEKP  (J3,J4) 

00001990 

ID2  *  ID (2)  /  6 

00002000 

ENDIF 


IF  (CHEFS (J3,J4)  .  EQ.  0)  THEN  DO 
JSEG  =  ID3  =  ID { 3)  =  1 
ELSE  DO 

JSEG  =  6  ;  ID  (3)  =  6  *  CHEKS(J3,J4) 

ID3  =  ID  (3)  /  6 

ENDIF 

:2  WRITE  (3,2000)  ITEST  (ID  ( 1)  )  ,  ID(1),  ID  (2)  , 

:2  S  (  (PANEL(J3,J4,J5,J6)  ,  J6  =  1,  JPAN)  ,  J5  =  1,  ID2) 

'.2  WRITE  (3,2100)  ID  (3),  ( (SEGMET  ( J3 ,  J4,  J5 ,  J6)  , J6  =  1,  JSEG)  ,  J5=  1 ,  ID3) 
:3  "WRITE  (IDISF2)  ID,  ( (P  ANEL  ( J3,  J4 ,  J5,  J6)  ,  J6=  1 ,  JPAN)  ,  J5=  1  ,ID2)  , 

:3  *  (  (SSGHET  (J3,J4,  J5,J6)  ,J6=1,JSEG)  ,J5=1,ID3) 

1000  ‘  CONTINUE 

STOP 

>00  PRINT, 'THE  VALUE  OF  THE  FLAG  IS  L. T.  ZERO  - ERROR  - • 

STOP 

100  PORHAT  ( 1 H 1 ,  10X,  •FOPCES  PROH  ATB  HOBEL*/) 

101  PORHAT  (44X,A10,2X,I2,2X,A3,19X,I2,2X,A3) 

102  PORHAT  {/,  20X,  12,  40X,  A3,  15X,  12,  40X,  A3) 

103  PORHAT  (44X,  A10) 

1416  PORHAT  (180  (/)  ,A10) 

113  PORHAT  (///A 10) 

114  PORHAT  (11 (/) , A10) 

115  FORMAT  (F9.3,9X,2F9.2,9X,3F8.3,9X,  2F9.2,9X,3F8.3) 

’23  PORHAT  (1 0  (/)  ,A10) 

125  PORHAT  (F9. 3,9X,2P9. 2,9X, 3P8.3, 2X, 3F8.3) 

’000  PORHAT  (1H  ,  A4,  215,  5X,  F5.0,  5F12.3) 

2100  PORHAT  (1H  ,  114,  5X,  F5.0,  5F12.3) 

2200  PORHAT  (»-TIHE  =*,F6.3,'  HSEC«) 

END 

5DATA 

5END 

1STOP 

V*PW=BONE 

'/ 


00002010 

00002020 

00002030 

00002040 

00002050 

00002060 

00002070 

00002080 

00002090 

00002100 

00002110 

00002120 

00002130 

00002140 

00002150 

00002160 

00002170 

00002180 

00002190 

00002200 

00002210 

00002220 

00002230 

00002240 

00002250 

00002260 

00002270 

00002280 

00002290 

00002300 

00002310 

00002320 

00002330 

00002340 

00002350 

00002360 

00002370 


VALF  ,?A  JOB  nU.D08.AQ0221,BRANDEAU,H={2,0) 

'/  EXEC  WATFIV 

'/GO. FT06F00 1  DD  DSN=DU.D08. AQO 221 . BRAND  BAD . ATBDATA. ONE, DISP* SHE 
'/GO. FT07F001  DD  DSN=DU.D08. AQO 2 2 1 . BP  AND EAU . ATBDATA. TWO , DISP=SHR 
'/GO. FT09F001  DD  DSN=DU. D08. AQ0221 . BR AND EAU . ATBDATA . FINAL, DISP=SHR 
'/GO. SYSIN  DD  * 

;job 

:$OPTIONS  NOEXT, NOCHECK, CCOMP=0 

REAL  VE  H  ( 1  3)  ,  FORCE(8,30),  PANEL(8,24),  INSEG(8,24)f  TIME 
PEAL  DUMMY  (7)  ,  OUTSEG  (8,72)  ,  POSIT  (24)  ,  CNSTRN(8,3) 

INTEGER  NT,  LIMB ,  ID  (8,3)  ,  KT(8) 

READ  (6)  NT,  DLT,  KC0N4 ,  (POSIT  (I),  1=1,  KCON4) 

WRITE  (9)  NT,  DLT,  KCON4,  (POSIT(I),  1=1,  KCON4) 

:  READ  AND  WRITS  SEGMENT  *,  HASS,  INERTIA,  AND  SEMI-MAJOR  AXES. 

DO  10  I  =  1,  8 

HEAD  (6)  J,  DUMMY 
10  WRITE  (9)  J,  DUMMY 

BEGIN  LOOP  FOR  ALL  REHAINING  DATA  TO  BE  CCHBINED 

DO  100  I  =  1,  NT 
READ  (€)  YEH 
READ  (7)  TIME 
PRINT, *  TINS  =  »,TIHE 
WRITE  (9)  VEH 

READ  IN  ALL  DATA  FOR  THIS  TIHE  STEP 
DO  40  J  =  1,  8 

READ  (6)  LIMB,  K3,  (FORCE(J,K),  K  =  1 ,30)  ,  (CNSTRN  (J  ,K)  ,K=  1 ,  K3) 
KT  (J)  =  K3 

CONVERT  ROTATIONS  FROM  DEGREES  TO  RADIANS 
DO  20  K  =  7,9 

20  FORCE  (J ,  K)  =  FORCE  (J,K)  /  57.29578 

READ  (7)  ID  (J,  1)  ,  ID2,  ID3,  (PANEL  (J,  KK)  ,  KK  =  1,  ID2)  , 

S  (INSEG  ( J,  KK)  ,  KK  =  1,  ID3) 

ID  (J,2)  =  ID2  ;  ID  (J,  3)  =  ID3 
IF  (LIMB  .NE.  ID(J,1))  THEN  DO 

PRINT , ' LI  MB  NOS.  NOT  EQUAL  AT  STEP  *,J 
PRINT, 'LIMB  NO.  FROM  UNIT  6  =',LIMB 
PRINT, 'LIMB  NO.  FROM  UNIT  7  =  ',ID(J,1) 

PRINT, 'STOPPING  HOW* 

STOP 

ENDIF 

40  CONTINUE 


00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250 

00000260 

00000270 

00000280 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 


COMBINE  DATA  FOB  EACH  CONTACTED  SEGMENT  SO  THAT  RELATIVE 
VELOCITIES  CAN  BE  CALCULATED.  INCREASE  SEGMENT  DATA  FROM 
:  6  TO  18  ITEMS. 

DO  60  J  =  1,  8 

ID  3  *  ID  ( J,  3)  ;  ID2  =  ID(J,2);  OUTSEG(J,1)  =  0.0 
IF  (ID3  .  EQ.  1)  GO  TO  55 
DO  50  K  =  1,  ID3,  6 

NUM  =  IFIX(INSEG  (J,K)  ) 

K1  =  3  *  k  -  2  ;  K2  =  K1  *  17  ;  K3  =  K 1  +  6  ;  K4  =  Kl  ♦  5 
DO  44  L  =  FI,  FU 

44  OUTSEG  (J,L)  =  INSEG  (J,  K+L-K1) 

IF  (NUM  .SQ.  0)  THEN  DO 

DO  *5  L  *  K3,  K4  • 

45  OUTSEG  (J, L)  =0.0 
ELSE  DO 

:  FIND  PROPER  LOCATION  IN  SEGMENT  "NUM"  FILE 
ID2  =  1 

WHILE  (IFIX  (INSEG  (NUN,ID2)  )  .NE.  J)  DO 
ID2  =  ID 2  +  6 
ENDWHILE 

DO  48  L  =  1,  : 

48  OUTSEG  (J,KH-L+5)  =  FORCE  (  NUM  ,L+3) 

DO  49  L  =  15,  17 

49  OUTSEG  (J, Kl *L)  =  INSEG  (NUM,ID2+L-1  2) 

ENDIF 

50  CONTINUE 

ID3  =  ID3  /  6  *  18 

:  NEGATE  "OTHER  FORCES"  FOR  EQUILIBRIUM 
55  K3  =  KT  (J) 

WRITE  (9)  ID  (J,1)  ,  ID2,  ID3,  K3,  (FORCE  (J,KK)  ,  KK  =  1 ,  30)  , 
$  (PANEL  (J,KK)  ,KK  =  1,  ID2)  ,  (OUTSEG  (J, KK)  ,  KK  =  1,  ID3) 

$  ,  (-CNSTRN  (J,KK)  ,  KK  =  1,  K3) 

:2  WRITE  (3,  300)  ID  (J,  1)  ,  ID2,  (PANEL  (J, KK)  ,  KK  =  1,  ID2) 

:2  WRITE  (3,200)  ID3,  (OUTS  EG  (J,  KK)  ,  KK  =  1,  ID3) 

:2  I?  (K 3  .GT.  1)  WPITE  (3,150)  K3,  (-CNSTRN  (J,  KK)  ,  KK  =  1,K3) 

60  CONTINUE 
100  CONTINUE 
STOP 

:2300  PORMAT  ('  LIMB  =',  12,13, 5X,  (2  <F5.  0 , 5F8.  2 , 5X) )  ) 

:2200  PORMAT  (»  SEGS* ,13 , (F5.0, IX,  2F6.2,  2X,  15F7.1)) 

:2150  FORMAT  ('EXTRA  CONTACTS  14,  3F12.4) 

END 

:DATA 

IBRD 


00000510 

00000520 

00000530 

00000540 

00000550 

00000560 

00000570 

00000580 

00000590 

00000600 

00000610 

00000620 

00000630 

00000640 

00000650 

00000660 

00000670 

00000680 

00000690 

00000700 

00000710 

00000720 

00000730 

00000740 

00000750 

00000760 

00000770 

00000780 

00000790 

00000800 

00000810 

00000820 

00000830 

00000840 

00000850 

00000860 

00000870 

00000880 

00000890 

00000900 

00000910 

00000920 

00000930 

00000940 

00000950 

00000960 

00000970 

00000980 

00000990 

00001000 


:$CPTIOSS  DEC,CCOMP=1 
:  PROGRAM  NAME  :  CONTACT 

PEAL  VSH(12),  FORCE  (30)  ,  PANEL(24),  SEGMET  (72)  ,  PLANE(33) 

REAL  ABC  ( 24)  ,  XY2{24),  TIME,  INERT  (24)  ,  WEIGHT  (8)  ,  D1  (3) 

REAL  RADS  (3),  OMEGA  (3),  THOLD  (3)  ,  OMEGA2  (3)  ,  RADS2(3) 

REAL  XROLD  ( 3)  ,  XTSHP(3),  TEMP(3),  WORK  (6)  ,  POSIT  (24),  CNSTRN  (6) 
CHARACTER  SEG*3(8)  /• RUL* , »RLL» , » LUL* , • ILL* , * RU A’ , » RLA' , • LUA • , 

*  ' LLA  *  / 

INTEGER  NT,  LIMB,  NPAN,  NS  EG ,  KT,  JS,  JS  2,  JS3,  LIMBO,  LIMBI,  KKK 
INTEGER  KK,  J1,  J2,  J,  I,  K,  ICHEK,  KS,  KP,  NUM,  NUMQ,  NCON,  KCON4 
COMMON  /TRANS/  ICHEK 
INTEGER  YES  (8) 

DATA  YES  /8  *  0  / 

-REAL  GINCH  /  386.0888  / 

DATA  PLANE  /-. 1 1 04 , 0. 0 ,-. 9 939,  . 8 744 , 0. 0, 2249, 

*  0. 0,0.0, -1.0,  .9745, 0.0, -.2245,  . 919 1 , 0. 0 ,. 3939 , 

*  -.0499, 0.0, .9988,  -  1 . 0, 0. 0, 0. 0,  0. 0, -  1 . 0 , 0. 0 , 

*  0.0, 1.0, 0.0,  .9720,0.0,-. 2350,  -. 6428, 0. 0,-. 7660  / 

DATA  D 1  /8. 8,  7.93,  8.8,  7.93,  5.44,  7.94,  5.44,  7.94  / 

EQUIVALENCE  (RADS,  F0RC£(7)},  (OMEGA, FORCE  (10) ) 

READ  (5)  NT,  DLT,  KCON4,  (POSIT(I),  1=1,  KCON4) 

IP  (KCON4  .EQ.  1)  GO  TO  4 

ROTATE  "OTHER  CONSTRAINT  FORCE"  LOCATION  TO  X-AXIAL  COORDINATES 
WITH  THE  ORIGIN  AT  THE.  PROXIMAL  JOINT 

DO  3  I  a  1,  KCON4 ,  4 

CALL  CHANGE  (POSIT,  1+1,  24,  D1 ( POSIT  (I) ) ) 

3  YES (POSIT (I))  =1+1 

4  WRITE  (9)  NT,  DLT,  KCON4,  (POSIT (I),  1=1,  KCON4) 

WRITE  (9)  D1(1),  D1  (2)  ,  D1  (5)  ,  D1(6) 

READ  SEGMENT  WEIGHTS,  INERTIAS,  AND  SEMI-MAJOR  AXES  IN  LOCAL  X-AXIAL 
SEGMENT  COORDINATES. 

J1  =  -2 
DO  5  J  =  1,  8 
J1  =  J1  +  3 
J2  =  J 1  ♦  2 

HEAD  (5)  I, WEIGHT  (J),  (INERT  (K),  K=J1,J2),  (ABC  (K)  ,  K  =  Jl,  J2) 
WRITE  (9)  I,  WEIGHT  (J)  ,  (ABC(K),  K=J1,J2) 

SQUARE  SEMI-MAJOR  AXES  OP  ELLIPSOIDS  POR  LATER  USE 

DO  5  K  =  Jl,  J2 

5  ABC  (K)  =  ABC  (K)  *  ABC  (K) 

DO  350  KKK  =  1,  NT 
READ  (5)  TIME,  7EH 


00000010 

00000020 

00000030 

00000040 

00000050 

00000060 

00000070 

00000080 

00000090 

00000100 

00000110 

00000120 

00000130 

00000140 

00000150 

00000160 

00000170 

00000180 

00000190 

00000200 

00000210 

00000220 

00000230 

00000240 

00000250 

00000260 

00000270 

00000280 

00000290 

00000300 

00000310 

00000320 

00000330 

00000340 

00000350 

00000360 

00000370 

00000380 

00000390 

00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000460 

00000470 

00000480 

00000490 

00000500 


U)  Ui  UJ  U  U  I  |  u  U  to  Ul  u 


1  PRINT,  ’TIME  =»,TIME 

WRITS  (9)  Tins 

DO  350  KK  =  1,  8 

READ  (5)  LIMB,  NPAN,  NSEG,  NCON,  FORCE,  (PANEL  (I)  ,  1=1,  NPAN) 
*  ,  (SEGMET  (I)  ,  I  =  1 ,  NSEG)  ,  (CNSTRN(I),  1=1,  NCON) 

LIMBO  =  3  *  LIMB  -  3 

3  PRINT,  *  **+**************************************■********  • 

3  PRINT, ’SEGMENT  IS  #’,LIMB,’  OR  ' ,S EG  (LIMB) 

3  PRINT, **************************************************’ 

3  PRINT, ’ROTATION  ANGLES  (RADIANS)’ 

3  WRITE  (3,500)  (?ORCF(I),I  =  7,  9) 

ICHEK  =  -1 
•KS  =  0 

KP  =  0 

IP  (NSEG  .EQ.  1)  GO  TO  50 
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SEGMENT  -  SEGMENT  CONTACT 

uu  u  uuo  /  u 

00000680 
.AAA  A  At Q  A 

SEGMENT  CONTACT  POINTS  IN  SEGMENT  COORDINATES 

00000700 

00000710 

KS  =  NSEG  /  18 

00000720 

DO  40  J  =  1,  KS 

00000730 

JXYZ  =  6  *  J  -  5 

00000740 

JS  =  1  ♦  (J-1)  *  18 

00000750 

JS2  =  JS  ♦  1 

00000760 

JS3  =  JS  ♦  2 

00000770 

NTTM  =  I FIX  (SEGMET  (JS)  ) 

0000078C 

3  PRINT,’  * 

00000790 

3  PRINT, ’SEGMENT’ , LIMB, •  CONTACTING  SEGMENT’ 

,  NU  M 

00000800 

3  PRINT, ' NORMAL  PORCE  = ’ ,SEG  MET  (JS2) 

>  00000810 

3  PRINT, ’FRICTION  FORCE  =* ,  SEGMET (JS3) 

00000820 

3  PRINT, ’CONTACT  AT  ATB  LOCAL  XYZ (SEGMET (JS  +  I) ,  I  =-  3, 

5) 

00000830 

- - - 

00000840 

NORMAL  CONTACT  FORCES 

00000850 

— 

00000860 

00000870 

FIND  NORMAL  TO  ELLIiSOID  AT  POINT  OF  CONTACT 

00000880 

00000890 

BOLD  =0.0 

00000900 

DO  10  I  =  1,  3 

00000910 

HOLD2  =  SEGMET (JS 3+1)  /  ABC(LIMB0+I) 

00000920 

TEMP (I)  =  HOLD2 

00000930 

WORK (I*3)  =  HOLD 2 

00000940 

10  HOLD  =  HOLD  ♦  HOLD 2  *  HOLD 2 

00000950 

HOLD  =  SQRT  (HOLD) 

00000960 

00000970 

MULTIPLY  UNIT  NORMAL  BY  NORMAL  FORCE  SCALAR  - 

STORE  IN  "WORK  (1—3)  " 

00000980 

00000990 

DO  15  I  =  1,  3 

00001000 

15  WORK  (II  =  -TEMP(I)  *  SEGMET  (JS2)  /  HOLD 
IF  (SEGMET  (JS3)  .E Q.  0.0)  GO  TO  39 

FEICTION  CONTACT  FORCES 


:  ASSIGN  DATA  FOR  SEGMENT  "NON"  TO  STORAGE  VECTORS 

DO  20  I  =  1,  3 
JSI  »  JS  ♦  I 

XTEMP  (I)  =  SEGMET  (JSI  +  14} 

THOLD  (I)  =  SEGMET  (JSI  +  5) 

RADS2  (I)  =  SEGMET  (JSI  ♦  8) 

•  OMBG  A2  (I)  =  SEGMET  (JSI  +  11) 

20  TEMP  (I)  =  SEGMET  (JSI  ♦  2) 

:  CONSTROCT  RELATIVE  VELOCITY  VECTOR  IN  SEGMENT  "LIMB"  COORDINATES 
:  FOR  SEGMENT  "LIMB"  CONTACTING  SEGMENT  "NOM". 

:  VECTORS  "OMEGA2",  "P.ADS2",  "THOLD"  &  "ITEM!"  CONTAIN  DATA  FOR 
:  SEGMENT  "NOM". 

CALL  CROSS  (OMEGA2,  XTEMP,  X HOLD) 

:9  PRINT, 'W2  CROSS  R2=', XHOLD 

CALL  ROT  (XHOLD,  1,  3,  HADS2,  -1) 

:9  PRINT, 'SAME  VECTOR  ROTATED  INTO  INERTI AL • , XHO LD 
ICHEK  =  -1 

CALL  CROSS  (OMEGA,  TEMP,  XTEMP) 

:9  PRINT, 'W  CROSS  R=’.,  XTEMP 

DO  25  I  =  1,  3 

25  THOLD (I)  =  FORCE (I >3)  -  THOLD (I)  -  XHOLD  (I) 

CALL  ROT  (THOLD,  1,  3,  FADS,  +1) 

DO  30  Z  *  1V  3 

THOLD  (I)  *  THOLD  (I)  ♦  XTEMP (I) 

30  TEMP  (I)  =  WORK  (1  +  3) 

:9  PRINT, 'RELATIVE  ‘ELOCITY  VECTOR  IN  LOCAL’, THOLD 

:  PROJECT  VELOCITY  VECTOR  ONTO  TANGENT  PLANE  TO  SEGMENT  "LIMB" 

CALL  VECTOR  (TEMP,  THOLD,  XHOLD) 

:9  PRINT, 'SAME  PROJECTED  INTO  PLANE  OF  SEG.  », XHOLD 

:  HOLTIPL Y  PROJECTED  VELOCITY  BY  FRICTION  SCALAR  AND  COMBINE 
:  ALL  FORCES  INTO  "WORK  (1-3)" 

DO  35  I  =  1,  3 

35  WORK  (I)  *  WORK  (I)  +  XHOLD  (I)  ♦  SEG  BET  (JS  3) 

I  CONVERT  SEGMENT  Z-AXIAL  COORDINATES  TO  X-AXIAL 

39  XTZ(JXYZ)  *  SEGMET  (JS  +  5)  ♦  DI(LIMB) 
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XYZ  ( JXYZ  +  1)  =  SHGMET  (JS+3) 
XYZ  ( JXYZ  +  2)  =  SEGMET  (JS+4) 
XYZ  (JTYZ  +  ?)  =  WORK  (3) 

XYZ  (JXYZ  +  4)  =  WORK  (1) 

XYZ  ( JXYZ  +  5)  =  WORK  (2) 

40  CONTINUE 

50  IF  (NPAN  .EQ.  1)  GO  TO  130 


SEGHENT  -  PANEL  CONTACT 


:  PANEL  CONTACT  DATA  IS  IN  INERTIAL  COORDINATES. 

XP  =  NP AN  /  6 
■  K 1  =  KS  ♦  1 
K2  -  K 1  +  KP  -  1 
DO  120  J1  =  K 1 ,  K2 
JXYZ  =  6  *  J1  -  5 
J  =  J1  -  KS 
JS  =  1  ♦  (J-1)  *  6 
JS2  =  JS  ♦  1 
JS3  =  JS  ♦  -2 

;  CORRECT  PANEL  CONTACT  POINT  FOR  VEHICLE  NOTION 

DO  55  12  *  1,  KP 
DO  55  I  *  1,  3 
NON  »  6  *  (J-1)  ♦  I  +  3 

55  PANEL(NUK)  =  PANEL  (NUN)  ♦  VEH(I) 

:3  NUN  =  IFIX  (PANEL  (JS)  ) 

:3  PRINT, '  • 

:3  PRINT, ' SEGMENT ' , LI  MB, '  CONTACTING  PANEL’, NON 

:3  PRINT, ’NORMAL  FORCE  =' , PAN  EL ( JS 2) 

:3  PRINT, ’FRICTION  FORCE  *’ , PANEL (JS3) 

:3  PRINT, ’CONTACT  AT  INERTIAL  XYZ' ,  (PANEL  (JS  +  I)  ,  I  =  3,  5) 

:  FIND  VECTOR  PRON  CENTER  OF  ELLIPSOID  TO  CONTACT  POINT 
I  IN  INERTIAL  COORDINATES. 

DO  6  0  I  =  1 ,  3 

TEMP (I)  *  PANEL (JS3+I)  -  FORCE (I) 

60  THOLD(I)  *  TEMP (I) 

:5  PRINT, 'INERTIAL  CONTACT  VECTOR  BEFORE  TRANSFORMATION* 

:5  WRITE  (3,500)  TEMP 

:  TRANSFORM  CONTACT  VECTOR  FROM  INERTIAL  TO  SEGMENT  COORDINATES 

CALL  ROT  (TEMP,  1,  3,  RADS,  1) 

:5  PRINT, 'AFTER  TRANSFORMATION  TO  SEGHENT  LOCAL' 

:5  WRITE  (3,500)  TEMP 

:5  PRINT, ’  • 
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V£>  SO  VO  VJO  <X> 


STORE  CONTACT  POINT  IN  X-  AXI  A  L  COORDINATES 

XYZ(JXYZ)  =  TEMP  (3)  +  D1  ("LIME) 

XYZ  ( JX  Y  Z  + 1 )  =  TBMP  (1) 

XYZ (JXYZ+2)  =  TBMP (2) 

NORMAL  CONTACT  FORCES 

NUM  =  IFIX  (PANEL  (JS) ) 

NUMO  =  3  *  NOM  -  3 
DO  70  I  =  1,  3 

TEMP  (I)  =  PLANE  (  NUMO+I) 

•  WORK  (1  +  3)  =  0.0 

XTEMP  (I)  =  OMBGA (I) 

70  WORK  (I)  =  TEMP  (I)  *  PANEL(JS2) 

IF  (PANEL  (JS3)  .EQ.  0.0)  GO  TO  95 

FRICTION  CONTACT  FORCES 


FIND  VELOCITY  VECTOR  OF  CONTACT  POINT  IN  INERTIAL  COORDINATES. 

CALL  ROT2  (XTEMP,  1,  3,  -1) 

CALL  CROSS  (XTEMP,  THOLD,  XHOLD) 

PRINT, 'ANGULAR  VELOCITY  OF  C.G. '  ' 

WRITE  (3,500)  XTEMP 

PRINT, 'LINEAR  VELOCITY  DUE  TO  ROTATION  —  OMEGA  CROSS  R' 

WRITE  (3,500)  XHOLD 

PRINT, 'TRANSLATIONAL  VELOCITY' 

WRITE  (3,500)  (FORCE  (I) ,  I  =  4,  6) 

DO  80  I  =  1,  3 

80  XHOLD (I)  =  XHOLD (I)  +  FORCE  (1+3) 

6  PRINT, 'TOTAL  VELOCITY  VECTOR* 

6  WRITE  (3,500)  XHOLD 

PROJECT  VELOCITY  VECTOR  ONTO  TANGENT  PLANE 

CALL  VECTOR  (TEMP,  XHOLD,  THOLD) 

6  PRINT, 'UNIT  VECTOR  OF  VELOCITY  PROJECTED  ONTO  TANGENT  PLANE* 

6  WRITE  (3,500)  THOLD 

MULTIPLY  NEG.  UNIT  VECTOR  BY  FRICTION  SCALAR  TO  GET  FRICTION  VECTOR 
DO  90  I  =  1,  3 

90  WORK  (1  +  3)  =  -THOLD  (I)  *  PANEL  (JS3) 

CREATE  TOTAL  FORCE  VECTOP  IN  INERTIAL  COORDINATES 
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in  in  id 


95  DO  100  1=1,3 
100  TEMP  (I)  =  WORK  (I)  +  WORK  {1+3) 

!6  PRINT, ' TOTAL  FORCE  VECTOR  IN  INERTIAL  COORDINATES* 

!6  WRITE  (3,500)  TEMP 

:  TRANSFORM  TOTAL  FORCE  VECTOR  TO  SEGMENT  Z-AXIAL  COORDINATES. 
CALL  ROT2  (TEMP,  1,  3,  1) 

PRINT, 'AFTER  CHANGING  TO  X-AXIAL  COORDINATES’ 

WRITE  (3,500)  TEMP 
PRINT, •  • 

CONVERT  SEGMENT  Z-AXIAL  COORDINATES  TO  X-AXIAL. 

XYZ  (JXYZ+3)  =  TEMP (3) 

XYZ  (JXTZ  +  4)  »  TEMP  (1) 

XYZ  (  JXYZ  +  5)  =  TEMP  (2) 

120  CONTINUE 


JOINT  FORCES  S  TORQUES;  LINEAR  6  ANGULAR  ACCELERATIONS 


:  IN  BOTH  INERTIAL  AND  LOCAL  COORDINATES 
;  LOCAL  ;  ANGULAR  VELOCITIES,  ACCELERATIONS,  U2  ARRAY 
:  INERTIAL  ;  JOINT  FORCES  S  TORQUES,  LINEAR  ACCEL.,  U1  ARRAY 

130  KT  =  KS  +  KP 
KT6  »  KT  *  6 

:  TRANSFORM  INERTIAL  VECTORS  TO  SEGMENT  X-AXIAL  LOCAL  VECTORS. 

:4  PRINT, 'W, FORCES,  TORQUES,  ALPHA,  ACCEL.,  U1,  U2  IN  ORIG.  COORDS’ 
J4  WRITE  (3,600)  (FORCE  (I)  ,1=  10 ,30) 

DO  134  J  =  13,  25,  3 

IF  (J  .EQ.  19)  GO  TO  134 

CALL  ROT  (FOTRCE,  J,  30,  RADS,  +1) 

CALL  CHANGE  (FORCE,  J,  3  0,  0.0) 

134  CONTINUE 

:  ROTATE  LOCAL  VECTORS  TO  X-AXIAL 

DO  136  J  =  10,  28,  9 
136  CALL  CHANGE  (FORCE,  J,  30,  0.0) 

:4  PRINT, 'SAME  ITEMS  IN  SEGMENT  X-AXIAL  COORDINATES' 

:4  WRITE  (3,600)  (FORCE(J),  J  =  10,  30) 

IF  (YES  (LIMB))  140,  140,  138 

:  TRANSFORM  "OTHER  CONSTRAINT  FORCE"  TO  SEGMENT  LOCAL  X-AXIAL 
Z  1 

138  CALL  ROT  (CNSTRN ,  1,  6,  RADS,  +1) 

CALL  CHANGE  (CNSTRN,  1,  6,  0.0) 
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on  cc  cc.  tr 


laO  KS  =  YES  (LI SB) 

IF  (KT  .LT.  1  .AND.  KS  .LT.  1)  GO  TO  150 

CORRECT  01  £  02  ARRAYS  TO  ALLOW  FOR  POINT  FORCE  APPLICATIONS. 

DO  142  I  =  1,  3 
142  XHOLD(I)  =  0.0 

IF  (KS  .LT.  1)  GO  TO  144 

SOM  TOHQOES  ABOOT  C.G.  DOE  TO  OTHER  CONSTRAINT  PORCES. 

XHOLD { 1 )  =  POSIT  (KS  +  1)  *  CNSTRN(3)  -  POSIT  (KS  +  2)  *  CNSTRN(2) 
XHOLD(2)  =  - (POSIT (KS) -D1 (LIMB)) *CNSTRN  (3)  ♦  POSIT  (KS+2) *CNSTRN  (1) 

•  XHOLD ( 3)  =  (POSIT  (KS)  -Dl  (LIMB) }  *CNST8N  (2)  -  POSIT  (KS  +  1)  *CNSTRN  ( 1) 
IF  (KT  .LE.  1)  GO  TO  146 

SOM  TORQOES  ABOOT  C.G.  DOE  TO  CONTACT  FORCES. 

144  GMASS  =  GINCH  /  WEIGHT  (LIMB) 

DO  145  J  =  1,  KT 

I  ■  6  *  J  -  5 

XHOLD  ( 1)  =  XHOLD  (1 )  +  XYZ  (1+1)  *  XYZ  (1  +  5) 

S  -  XTZ  (1*2)  *  XYZ  (1+4) 

XHOLD ( 2)  =  XHOLD  (2)  -  (XYZ  (I)-D1  (LIMB) )  *  XYZ  (I ♦5) 

$  +  XYZ  (1  +  2)  *  XYZ  (1+3) 

XHOLD ( 3)  =  XHOLD  (3)  +  (XYZ  (I)-D1  (LIMB) )  *  XYZ  (1+4) 

*  -  XYZ  (1  +  1)  *  XYZ  (1+3) 

DO  145  JJ  ■  3,  5 

FORCE (JJ  +  22)  =  FORCE  (JJ+22)  -  GMASS  *  XYZ(I+JJ) 

145  CONTINOE 

146  CONTINOE 

PRINT, '02  ARRAY  IN  X-AXIAL  COORDINATES* 

WRITE  (3,500)  (FORCE ( J) ,  J  ■  28,30) 

PRINT, » TORQOES  DOE  TO  ALL  POINT  SEGMENT  FORCES  ABOOT  C.  G. • 

WRITE (3 ,500)  XHOLD 

ADD  T/I  TO  02  TO  PRODOCE  ADJOSTED  ACCELERATION. 

DO  148  J  s  1,  3 

148  FORCE  (J+27)  =  FORCE  (J+27)  +  XHOLD(J)  /  I NERT  ( J+LI51B0) 

8  PRINT, ' ADJOSTED  02* 

8  WRITE  (3,500)  (FORCE  (J)  ,  J  =  28,30) 

COMBINE  A  AND  01  INTO  LOCATION  OF  A. 

150  DO  152  J  *  22,  24 

152  FORCE  (J)  =  FORCE  (J)  -  FORCE  (J+3)  /  GINCH 
ORDER  CONTACTS  BY  INCREASING  X  VALOE  OF  CONTACTING  POINT 
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IF  (KT  .LB.  1)  GO  TO  200 
K  V  =  KT  -1 
DO  170  1=1,  KT1 
J1  =  f-  *  I  -  5 
HOLD  =  XYZ(Jl) 

II  =  I  ♦  1 
HUH  =  J1 

DO  153  J  =  II,  KT 
J2  »  6  *  J  -  5 

IF  (XYZ  (J2)  .GE.  HOLD)  GO  TO  158 

NON  =  J2 

HOLD  =  XTZ  (J2) 

158  CONTINUE 

•  IF  {NOH  .EQ.  J1)  GO  TO  170 
J1  «  J1  -  1 
NO  PI  =  NOP!  -  1 
DO  16  0  K  =  1,  6 
J3  -  6  *  K  -  5 
HOLD  =  XYZ  (J  1  ♦  K) 

XYZ(J1+K)  =  XYZ  (NUH+K) 

160  XYZ ( NUM+K)  =  HOLD 

170  CONTINUE 

200  IF  (KT)  205,  205,  210 
205  KT  *  1 

XYZ  ( 1)  =0.0 
GO  TO  300 
210  KT  =  KT  *  6 
300  CONTINOE 

4  PRINT, ’CONTACT  FORCES' 

4  PRINT, 'NUMBER  OF  DATA  I TENS  =»,KT 
4  K2  =  KT  /  6 

4  IF  (KT  .GT.  1)  WRITE  (3,600)  (XYZ(I),  1=1,  KT) 

OOTPOT  RSSOLTS  TO  DISK 

WRITE  (9)  LIH3,KT,NCON,  (FORCE  (I)  ,1  =  10,24)  ,  (FORCE  (I)  , I  =  28,30), 
$  (XYZ (I)  ,1  =  1,  KT)  ,  (CNSTRN (I)  ,  1  =  1,  NCON) 

350  CONTINOE 
STOP 

500  FORMAT (1H  ,3F10.4) 

600  FORMAT  ( 1H  ,3  F10.4,' - *,3F10.4) 

END 

SUBROUTINE  POT  (V,  IZ,  N,  P,  K) 

SUBROUTINE  TO  TRANSFORM  VECTOR  V  THROUGH  P  RADIANS  INTO  ANOTHER 
COORDINATE  SYSTEM.  K  DETERMINES  WHETHER  THE  TRANSFORMATION  IS  FROM 
INERTIAL  TO  SEGMENT  OR  VICE-VERSA  : 

K  =  1  INERTIAL  TO  SEGMENT  TRANSFORMATION 

K  *-1  SEGMENT  TO  INERTIAL  " 
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REAL  5  (3,3)  ,  V  (»T)  ,  7N(3)  ,  P(3) 
COMMON  /TP.  A  NS/  ICHEK 


IF  (ICHEK 
ICHEK  =  1 
Cl  *  COS 


GO  TO  10 


COS  (P  (1)  ) 
COS  (P  (2) ) 
COS  (P  (3) ) 
SIN  (P  (1) ) 
SIN  (P  (2)  ) 
SIN  (P  (3)  ) 


IF  (K)  40,  60,  20 
20  DO  30  I  =  1,  3 
7N  (I)  =0.0 
DO  30  J  *  1,  3 

30  7N  (I)  *  VN  (I)  ♦  R  (I  ,  J)  *  V  (IZ1+J) 

GO  TO  60 

40  DO  50  I  =  1,  3- 
VN  (I)  =  0.0 
DO  50  J  =  1,  3 

50  VN  (I 1  =  VN{I)  ♦  R  (J,I)  *  V(J) 

60  DO  70  I  =  1,  3 
70  V(IZHI)  =  VN  (I) 

RETURN 

END 

SUBROUTINE  VECTOR  (  N,  V,  C) 

SUBROUTINE  TO  FIND  THE  VECTOR  PROJECTION  OF  VECTOR  V  ONTO  THE 
PLANE  WHICH  N  IS  NORMAL  TO.  THE  CROSS  PRODUCT  IS  TAKEN  TWICE 
AND  THE  RESULTANT  VECTOR  IS  NORMALIZED  TO  A  UNIT  VECTOR. 

BEAL  N  (3)  ,  V  (3)  ,  C(3),  D(3),  HOLD 

HOLD  =0.0 

CALL  CROSS  (N,  V,  C) 

CALL  CROSS  (C,  N,  D) 

DO  10  J  ■  1,  3 
10  HOLD  »  HOLD  ♦  D  (J)  *  D  (J) 

HOLD  «  SORT  (HOLD) 

IF  (HOLD  .EQ.  0.0)  RETURN 
DO  20  J  »  1,  3 


00004010 

00004020 

00004030 

00004040 

00004050 

00004060 

00004070 

00004080 

00004090 

00004100 


R  (1,1) 

= 

C2  *  Cl 

00004110 

R(2,  1) 

3 

S3  *  S2 

* 

Cl 

-  SI 

* 

C3 

00004120 

R  (3,1) 

3 

C3  *  S2 

* 

Cl 

♦  SI 

* 

S3 

00004130 

R  ( 1 , 2) 

3 

C2  *  SI 

00004140 

R  (2, 2) 

3 

S3  *  S2 

* 

SI 

♦  Cl 

* 

C3 

00004150 

R  (3, 2) 

3 

C3  *  S2 

* 

SI 

-  Cl 

* 

S3 

00004160 

P  (1,3) 

3 

-S2 

00004170 

R  (2,3) 

3 

S3  *  C2 

00004180 

R  (3,3) 

3 

C3  *  C2 

00004190 

ENTRY 

ROT 2  (V,  IZ 

t 

N, 

K) 

00004200 

IZ1  = 

IZ 

-  -1 

00004210 

00004220 

00004230 

00004240 

00004250 

00004260 

00004270 

00004280 

00004290 

00004300 

00004310 

00004320 

00004330 

00004340 

00004350 

00004360 

00004370 

00004380 

00004390 

00004400 

00004410 

00004420 

00004430 

00004440 

00004450 

00004460 

00004470 

00004480 

00004490 

00004500 


20  C<J)  =  D  ( J)  /  HOLD 
RETURN 
BSD 

SUBROUTINE  CROSS  (  A,  B,  C  ) 

SUBROUTINE  TO  TAKE  THE  CROSS  PRODUCT  OF  A  CROSS  B,  AND  RETORN 
IT  IN  VECTOR  C.  A  ,3  AND  C  ABE  ALL  VECTORS  OF  LENGTH  3. 

REAL  A  ( 3)  ,  Q  (3)  ,  C(3) 

C  ( 1)  =  A  (2)  *  B  (3)  -  A  (3)  *  B  (2) 

C  (2)  =  -A[1)  *  B  (3)  +  A  (3)  *  B  (1) 

C{3)  =  A  ( 1)  *  B  (2)  -  A  (2)  *  B(1) 

RETURN 
'  END 

SUBROOTINE  CHANGE  (TEMP,  I ,  N,  D) 

CHANGES  Z-AXIAL  TO  X-AXIAL  COORDINATES.  SHIFTS  NEW  X  BY 
D  TO  ALLOW  TRANSFER  OF  ORIGIN  TO  PROXIMAL  JOINT. 

REAL  TFSP(N) 
hold:  •-  TEMP  (I) 

HOLD2  =  TEMP  (I*1) 

HOLD 3  =  TEMP  (1  +  2) 

TEMP  (I)  =  HOLD3  ♦  D 
TEMP  (1+1)  =  HOLD1 
TEMP  (1  +  2)  *  HOLD2 
RETURN 
END 

SUBROUTINE  FOR  (K) 

DO  10  J  =  1,  K 
10  READ  (5) 

RETURN 

END 

SUBROUTINE  BACK  (K) 

DO  10  J  =  1,  K 
10  BACKSPACE  5 
RETURN 
END 


00004510 
00004520 
00004530 
00004540 
00004550 
00004560 
00004570 
00004580 
00004590 
00004600 
00004610 
00004620 
00004630 
00004640 
00004650 
00004660 
000  04670 
00004680 
00004690 
00004700 
00004710 
00004720 
00004730 
00004740 
00004750 
00004760 
00004770 
00004780 
00004790 
00004800 
00004810 
00004820 
00004830 
00004840 
00004850 
00004860 
00004870 
00004880 


FINDS  AREA  OT  30NF  AT  VARIOUS  LENGTHS 

00000010 

INTEGER  TITLE  (IS) 

00000020 

REAL  A ( 40)  ,  D  (40)  ,  DNCR3(40),  X(40),  1(40),  Z  (40)  ,  U, 

iP 

00000030 

READ  (1,100)  TITLE 

00000040 

WRITE  (3,150)  TITLE 

00000050 

PRINT, 'NORN  LENGTH,  LENGTH,  AREA* 

00000060 

PPINT, *  * 

00000070 

K  =  1 

00000080 

10  READ  ( 1  END=20)  D  (K)  ,  A  (K) 

00000090 

K  =  K  ♦  1 

00000100 

GO  TO  10 

00000110 

20  K  =  K  -  1 

00000120 

DO  30  I  =  1,  K 

00000130 

■  DNORH(I)  =  D(I)  /  D  (K) 

00000140 

30  WRITE  (3,200)  DNOPN  (I)  ,  D  (I)  ,  A  (I) 

00000150 

CALL  SPLINE  (K,  DNOPN,  A,  X,  Y,  Z) 

00000160 

00000170 

EVALUATE  3  11  POINTS  NORMALIZED 

00000180 

00000190 

PRINT, '  » 

00000200 

PRINT, •  « 

00000210 

PRINT, 'NORM  LENGTH  &  AREA  EVALUATED  BY  SEVAL' 

00000220 

PRINT, •  » 

00000230 

DO  40  I  »  1,  11 

00000240 

O  *  FLOAT  (1-1)  /  10.0 

00000250 

TEHP  =  SEVAL  (K,  O,  DNORN,  A,  X,  I,  Z) 

00000260 

40  RRITE  (3,200)  U,  TEHP 

00000270 

STOP 

00000280 

100  FORHAT  (18A4) 

00000290 

150  PORHAT  (*  ID  =  ',18A4) 

00000300 

200  FORHAT  (3F10.3) 

00000310 

END 

00000320 

SUBROUTINE  SPLINE  (N,  X,  Y,  E,  C,  D) 

00000330 

INTEGER  N 

00000340 

REAL  X  (N)  ,  Y  (N)  ,  B  (N)  ,  C(N),  D(N) 

00000350 

INTEGER  NH1,  IB,  I 

00000360 

REAL  T 

00000370 

NH1  *  N-1 

00000380 

IF  (  N  .LT.  2  )  RETORN 

00000390 

IF  {  N  ,LT.  3  )  GO  TO  50 

00000400 

D  (1)  =  1(2)  -  X  (1) 

00000410 

C  (2)  =  (Y  (2)  -  Y(1))/D(1) 

00000420 

DO  10  I  -  2,  NH1 

00000430 

D  (I)  =  X  (1*1)  -  X  (I) 

00000440 

B  (I)  »  2.*(D(I-1)  ♦  D  (I) ) 

00000450 

C(I*1)  »  (Y  (1  + 1 ) .  -  Y  (I)  )/D  (I) 

00000460 

C(I)  »  c  (1*1)  -  C(I) 

000004^0 

10  CONTINUE 

000004:0 

B  (1)  a  -D(1) 

00000490 

B  (N)  a  -D  (N-1) 

00000500 

C(1)  *  0. 

C(N)  »  0. 

IP  (  N  .  EQ.  3  )  GO  TO  15 

CO)  =  C  (3)  /  (X  (4)  -X  (2)  )  -  C  (2)  /  (X  (  3)  -X  (1) ) 

Cff»)  =  C(N-1)/{X{N) -X(N-2)  )  -  C{N-2)/{X(S-1)-X(M-3}) 

C(1)  =  C(1)  *D(1)**2/(X(4)-X(1)  ) 

C  (N)  =  -C  (N)  *D  (N-1)  **2/  (X  (N)  -X  (N-3) ) 

15  DO  20  I  =  2,  N 

T  *  D(I-1)/B(I-1) 

B  (I)  =  B(I)  -  T*D{I-1) 

C(I)  =  Cfl)  -  T*C  (I  —  1  > 

20  CONTINUE 

C  (M)  =  C(N)/B(N) 

DO  30  IB  ■  1,  N!11 
I  *  N-IB 

C  (I)  =  (Cfl)  -  D  (I)  *C  (I  +1 )  >  /B  (I) 

30  CONTINUE 

B  (N)  =  (Y{N)  -  Y  (Nfll)  J/DfNNI)  ♦  D (NB1) * (C  (SHI)  ♦  2.*C{N)) 
DO  40  I  =  1,  !»*1 

BCI)  =  (Y  (1  +  1)  -  Y  (I) )  /D  (I)  -  D  (I)  *  (C  fl+1)  +  2.*C(IJ) 
D(I)  =  fC{I*1)  -  C(I))/D(I) 

C  (I)  *  3.*C  (I) 

40  CONTINUE 

CfN)  =  3- *C  (N) 

D(N)  =  D  (N- 1) 

RETURN 

50  B  (11  *  {Y  (2)  -  Y  (1)  )  /  (X  (2)  -X  (1) ) 

C(1)  *  0. 

D  ( 1 )  =  0. 

B  (2)  *  B  ( 1) 

C (2)  *  0. 

D (2)  =  0. 

RETURN 

END 

BEAL  FUNCTION  SE7AL  (N,  0,  X,  Y,  B,  C,  D) 

INTEGER  N 

REAL  0,  X(N),  Y  (N)  ,  B  (N)  ,  C  (N)  ,  D  (») 

INTEGER  I,  0,  K 
PEAL  DX 
DATA  1/1/ 

IP  (  I  .GE.  N  )  I  =  1 
IF  (  0  .LT.  X(I)  )  GO  TO  10 
IF  (  0  .LE.  I  (1+1)  )  GO  TO  30 
10  I  =  1 
J  =  N+1 

20  PC  »  (I+J1/2 


00000510 

00000520 

00000530 

00000540 

00000550 

00000560 

00000570 

00000580 

00000590 

00000600 

00000610 

00000620 

00000630 

00000640 

00000650 

00000660 

00000670 

00000680 

00000690 

00000700 

00000710 

00000720 

00000730 

00000740 

00000750 

00000760 

00000770 

00000780 

00000790 

00000800 

00000810 

00000820 

00000830 

00000840 

00000850 

00000860 

00000870 

00000880 

00000890 

00000900 

00000910 

00000920 

00000930 

00000940 

00000950 

00000960 


IF 

( 

a 

.LT.  X  (K) 

)  o 

=  K 

00000970 

IF 

( 

u 

.GE.  X(X) 

)  I 

-  K 

00000980 

IF 

{ 

-i 

.GT.  1+1  1 

GO 

TO  20 

00000990 

30  DX 

St 

0 

-  X  (I) 

00001000 

\ 


Cl  0) 


;data 

•EMUF 

SEVAL  =  Y(I)  ♦  DX*  (B  (I)  ♦  DX*(C(I)  ♦  DX*D  (I) )  ) 

BETORN 

END 

#1  —  RIGHT  DISTAL  VS.  AVG.  INERTIA 

00001010 

00001020 

00001030 

00001040 

00001050 

'.0 

.226 

00001060 

0 

.092 

00001070 

l.  0 

.076 

00001080 

1.0 

.075 

00001090 

1.  0 

.074 

00001100 

5.0 

.069 

oooomo 

5.2 

.067 

00001120 

7.2 

.065 

00001130 

?.2 

•  .066 

00001140 

>.2 

.067 

00001150 

10.2 

.073 

00001160 

11.2 

.085 

00001170 

12.2 

.105 

00001180 

13.2 

.249 

00001190 

14.2 

1.818 

00001200 

DOUBLE  PRECISION  FUNCTION  ZEPOI N(AX,BX,F, TOL) 

DOUBLE  PRECISION  AX,3X,F,TOL 

A  ZERO  OF  THE  FUNCTION  F  (X)  IS  COMPUTED  IN  THE  INTERVAL  AX,BX  . 
INPUT.. 

AX  LEFT  ENDPOINT  OF  INITIAL  INTERVAL 

BX  RIGHT  ENDPOINT  OF  INITIAL  INTERVAL 

F  FUNCTION  SUBPROGRAM  WHICH  EVALUATES  F  (X)  POR  ANY  X  IN 

THE  INTERVAL  AX,BX 

TOL  DESIRED  LENGTH  OF  THE  INTERVAL  OF  UNCERTAINTY  OF  THE 
FINAL  RESULT  (  .GE.  0.0D0) 


OUTPUT.. 

ZEROIN  ABCISSA  APPROXIMATING  A  ZERO  OF  F  IN  THE  INTERVAL  AX,BX 


IT  IS  ASSUMED  THAT  F(AX)  AND  F(BX)  HAVE  OPPOSITE  SIGNS 
WITHOUT  A  CHECK.  ZEROIN  PETURNS  A  ZERO  X  IN  THE  GIVEN  INTERVAL 
AX,BX  TO  WITHIN  A  TOLERANCE  4*HACHEPS* ABS  (X)  +  TOL,  WHERE  HACHEPS 
IS  THE  RELATIVE  MACHINE  PRECISION. 

THIS  FUNCTION  SUBPROGRAM  IS  A  SLIGHTLY  MODIFIED  TRANSLATION  OF 
THE  ALGOL  60  PROCEDURE  ZERO  GIVEN  IN  RICHARD  BRENT,  ALGORITHMS  FOR 
MINIMIZATION  WITHOUT  DERIVATIVES,  PRENTICE  -  HALL,  INC.  (1973). 


DOUBLE  PRECISION  A , B, C,D , E, EPS , FA ,FB ,FC , TOL1 , XM ,P, Q,R , S 
DOUBLE  PRECISION  DABS,DSIGN 

COMPUTE  EPS,  THE  RELATIVE  MACHINE  PRECISION 

EPS  =  1.0D0 
10  EPS  =  EPS/2. 0D0 
TOL 1  =  1.0D0  ♦  EPS 
IF  (TOL  1  .GT.  1.0D0)  GO  TO  10 

INITIALIZATION 

A  =  AX 
B  =  BX 
FA  ■  F(A) 

FB  -  F  (B) 

BEGIN  STEP 


20  C  »  A 
FC  *  FA 


30  IF  (DABS(FC)  . G5.  DABS  (FB) )  GO  TO  40 
A  =  B 
B  -  C 
C  =  A 
FA  *  FB 
FB  =  FC 
FC  =  FA 

CONVERGENCE  TEST 

40  TOL 1  =  2. 0D0*EPS*D  ABS (B)  +  0.5D0*TOL 
XE  =  .  5*  (C  -  B) 

IF  (DABS  (XK)  .15.  TOLl)  GO  TO  90 
IF  (FB  .EQ.  0.0 DO)  GO  TO  90 

IS  BISECTION  NECESSARY 

IF  (DABS (E)  .LT.  TOLl)  GO  TO  70 
IF  (DABS  (FA)  .LE.  DABS  (FB)  )  GO  TO  70 

IS  QUADRATIC  INTERPOLATION  POSSIELE 

IF  (A  .NE.  C)  GO  TO  50 

LINEAR  INTERPOLATION 

S  =  FB/FA 
P  =  2.0D0*XN*S 
Q  =  I.ODO  -  S 
GO  TO  60 

INVERSE  QUADRATIC  INTERPOLATION 

50  Q  =  FA/FC 
R  a  FB/FC 
S  *  FB/FA 

P  *  S*  (2. 0D0*XH*Q*  (Q  -  R)  -  (B  -  A)*(H  -  I.ODO)  ) 
Q  =  (Q  -  1 . 0D0)  *  (R  -  1 . 0D0 )  *  (S  -  I.ODO) 

ADJUST  SIGNS 

60  IF  (P  .GT.  0.0D0)  Q  =  -Q 
P  =  DABS  (P) 

IS  INTERPOLATION  ACCEPTABLE 


IF  (  (2. 0D0*P)  .GE.  (3. 0D0*XM*Q  -  DABS  (TOLl *Q)  )  )  GO  TO  70 
IF  (P  .GE.  DABS (0. 5D0*E*Q) )  GO  TO  70 


E  =  D 
D  «  P/Q 
GO  TO  80 

BISECTION 

70  D  »  XM 
E  =  D 

COMPLETE  STEP 

80  A  =  B 
FA  =  FB 

‘IF  (DABS  (D)  .GT.  TOL1)  B  *  B  «■  D 
IF  (DABS(D)  .LE.  TOL1)  B  *  B  +  DSIGN(TOL1,  XM) 
PB  =  F (B) 

IF  (  (FB*(FC/DABS (FC) ) )  .GT.  0.0D0)  GO  TO  20 
GO  TO  30 

DONE 

90  ZEEOIN  =  B 
RETURN 
END 


SUBROUTINE  SPLINE  (N,  X,  Y ,  B,  C,  D) 

INTEGER  N 

DOUBLE  PRECISION  X(N),  Y  (N)  ,  B  ( N)  ,  C{N),  D  (N) 


THE  COEFFICIENTS  B  (I)  r  C  (I)  ,  AND  D(I),  1  =  1, 2,  ...,N  ARE  COMPUTED 
FOR  A  CUBIC  INTERPOLATING  SPLINE 

S(X)  =  Y  (I)  ♦  3  (I )  *  (X-X  (I)  )  +  C(I)*  (X-X(I) )  **2  ♦  D  (I)  *  (X-X  (I)  )  **3 

FOR  X(I)  ,LE.  X  .LE.  X(I+1) 

INPUT.. 

N  =  THE  NUMBER  OF  DATA  POINTS  OR  KNOTS  (N.GE.2) 

X  =  THE  ABSCISSAS  OF  THE  KNOTS  IN  STRICTLY  INCREASING  ORDER 

Y  =  THE  ORDINATES  OF  THE  KNOTS 

OUTPUT.. 

B,  C,  D  =  ARRAYS  OF  SPLINE  COEFFICIENTS  AS  DEFINED  ABOVE. 

USING  P  TO  DENOTE  DIFFERENTIATION, 

Y  (I)  =  S  (X  (I) ) 

B  (I)  =  SP(X(I)) 

C(I)  =  SPP(X(I))/2 

D  (I)  =  SPPP  (X  (I) )  /6  (DERIVATIVE  FRO!!  THE  RIGHT) 

THE  ACCOMPANYING  FUNCTION  SUBPROGRAH  SEVAL  CAN  BE  USED 
TO  EVALUATE  THE  SPLINE. 


INTEGER  NH1,  IB,  I 
DOUBLE  PRECISION  T 

NM1  =  N-1 

IF  (  N  .LT.  2  )  RETURN 
IF  (  N  .LT.  3  )  GO  TO  50 


SET  UP  TRIDIAGONAL  SYSTEM 

B  =  DIAGONAL,  D  =  OFFDI AGONAL ,  C  =  RIGHT  HAND  SIDE. 


D(1)  =  X  (2)  ~  X  (1) 

C  (2)  =  (1(2)  -  Y(1))/D(1) 

DO  10  1  =  2,  NM 1 

D(I)  =  X(I  +  1)  -  X  (I) 

B  (I)  «  2  . *  (D  (1-1)  +  D  (I) ) 

C  (1  +  1)  =  (Y  (1  +  1)  -  Y(I))/D(I) 
C(I)  =  C(IM)  -  C(I) 


10  CONTINUE 


END  CONDITIONS.  THIRD  DERIVATIVES  AT  X(1)  AND  X  (N) 
OBTAINED  FROM  DIVIDED  DIFFERENCES 

B  ( 1)  =  "D  ( 1 ) 

B  (N)  =  -D(N-1) 

C  ( 1)  =  0. 

C (N)  =  0. 

IF  (  N  ,EQ.  3  )  GO  TO  15 

C  (1)  =  C  (3)/(X(4)-X(2))  -  C  (2)  /  {X  {  3)  -X  (1)  ) 

C(N)  =  C(S-1)/(X(H)-X(N-2])  -  C  (N- 2)  /  (X  (N- 1)  -X  (N-3)  ) 

C(  1)  =  C(1)  *D(1)**2/(X(4)-X(1)) 

C(N)  =  -C  (N)  *D  (N-1)  **2/(X  (N)  -X  (N-3)  ) 

FORWARD  ELIMINATION 

15  DO  201  =2,  N 

T  =  D  (I-1)/B(I-1) 

B  (I)  =  B  (I)  -  T*D  (1-1) 

C(I)  =  C  (I)  -  T*C(I-1) 

20  CONTINUE 

BACK  SUBSTITUTION 

C  (N)  =  C(N)/B(N) 

DO  30  IB  =  1,  NM 1 
I  =  N— IB 

C  (I)  =  fC  Cl)  -  D(I)*C(I  +  1))/B(I) 

30  CONTINUE 

C(I)  IS  NOW  THE  SIGMA  (I)  OF  THE  TEXT 

COMPUTE  POLYNOMIAL  COEFFICIENTS 

B  (N)  =  (Y  (N)  -  Y  (NM  1) )  /D  (N  Ml)  «■  D  (NH1 )  *  (C  (NM1 )  +  2.*C(N)) 
DO  40  I  =  1,  NM1 

B  (I)  =  (Y  (1+1)  -  Y(I))/D(I)  -  D  (I)  *  (C  (1+1)  ♦  2.  *C  (I)  ) 

D  (I)  =  (C  (1  +  1)  -  C(I))/D(I) 

C(I)  =  3  .*C  (I) 

40  CONTINUE 

C  (N)  =  3  .  *C  (N) 

D  (N)  =  D  (N-  1) 

RETURN 

50  B  (1)  =  (Y  (2)  -  Y  (1) )  /  (X  (2)  -X  (1) ) 

C  (1)  =  0. 

D  ( 1)  =  0. 

B  (2)  *  B  (1) 

C (2)  =  0. 


D  (2)  = 
PETURN 
END 


DOUBLE  PRECISION  FUNCTION  SEVAL  (N,  U,  X,  Y,  B,  C,  D) 
INTEGER  N 

DOUBLE  PRECISION  U,  X  ((?}  ,  Y  (N)  ,  B  (N)  ,  C  (N)  ,  D  (N) 


THIS  SUBROUTINE  EVALUATES  THE  CUBIC  SPLINE  FUNCTION 

SEVAL  =  Y(I)  ♦  B  (I)  *  (TJ-X  (I)  )  ♦  C  (I)  *  {U-X  (I)  )  **2  *  D (I)  * (U-X (I) ) **3 

WHERE  X (I)  .LT.  0  .LT.  X(I+1) ,  USING  HORNER’S  ROLE 

IF  0  .LT.  X(1)  THEN  1=1  IS  USED. 

IF  U  .GE.  X  (N)  THEN  I  =  N  IS  USED. 

INPUT.. 


N  =  THE  NOHBER  OF  DATA  POINTS 

U  =  THE  ABSCISSA  AT  WHICH  THE  SPLINE  IS  TO  BE  EVALUATED 

X,Y  =  THE  ARRAYS  OF  DATA  ABSCISSAS  AND  ORDINATES 

B,C,D  =  ARRAYS  OF  SPLINE  COEFFICIENTS  COMPUTED  BY  SPLINE 

IF  U  IS  NOT  IN  THE  SAME  INTERVAL  AS  THE  PREVIOUS  CALL,  THEN  A 
BINARY  SEARCH  IS  PERFORMED  TO  DETERMINE  THE  PROPER  INTERVAL. 


INTEGER  I,  J,  K 
DOUBLE  PRECISION  DX 
DATA  1/1/ 


IF 

{  I  .GE. 

N  )  I 

=  1 

IF 

(  U  .LT. 

X(I1  ) 

GO  TO  10 

IF 

(  U  .LE. 

X(I+1) 

)  GO  TO  30 

BINARY 

SEARCH 

10  I  = 

1 

J  = 

N+1 

20  K  = 

(I  +  J)  /2 

IF 

(  U  .LT. 

X(K)  ) 

J  =  K 

IF 

(  U  .GE. 

X(K)  ) 

I  =  K 

IF 

(  J  .GT. 

1*1  ) 

GO  TO  20 

C  EVALUATE  SPLINE 


30  DX  =  U  -  X(I) 

SEVAL  »  Y(I)  ♦  DX*  (B  (I)  +  DX*(C(I)  ♦  DX*D(I))) 

RETURN 

END 


k 


non 


DOUBLE  PRECISION  FUNCTION  FMTN { AX, BX , F, TOL) 

DOUBLE  PRECISION  AX,3X,F,T0L 

AN  APPROXIMATION  X  TO  THE  POINT  WHERE  F  ATTAINS  A  MINIMUM  ON 
THE  INTERVAL  (AX,BX)  IS  DETERMINED. 


INPUT.. 

LEFT  ENDPOINT  OF  INITIAL  INTERVAL 
RIGHT  ENDPOINT  OF  INITIAL  INTERVAL 

FUNCTION  SUBPROGRAM  WHICH  EVALUATES  F(X)  FOR  ANY  X 
IN  THE  INTERVAL  (AX,BX) 

DESIRED  LENGTH  OF  THE  INTERVAL  OF  UNCERTAINTY  OF  THE  FINAL 
RESULT  (  .GE.  0.0D0) 


C  OUTPUT.. 

c 

C  FMIS  A3CISSA  APPROXIMATING  THE  POINT  WHERE  F  ATTAINS  A  MINIMUM 

C 

Z  THE  METHOD  USED  IS  A  COMBINATION  OF  GOLDEN  SECTION  SEARCH  AND 

C  SUCCESSIVE  PARABOLIC  INTERPOLATION.  CONVERGENCE  IS  NEVER  MUCH  SLOWER 
C  THAN  THAT  FOR  A  FIBONACCI  SEARCH.  IP  F  HAS  A  CONTINUOUS  SECOND 
C  DERIVATIVE  WHICH  IS  POSITIVE  AT  THE  MINIMUM  (WHICH  IS  NOT  AT  AX  OR 
C  BX),  THEN  CONVERGENCE  IS  SUPEBLINBAB,  AND  USUALLY  OF  THE  ORDER  OF 
C  ABOUT  1.324.... 

C  THE  FUNCTION  F  IS  NEVER  EVALUATED  AT  TWO  POINTS  CLOSER  TOGETHER 

C  THAN  EPS*ABS  (FHIN)  ♦  (TOL/3) ,  WHERE  EPS  IS  APPROXIMATELY  THE  SQUARE 
C  ROOT  OF  THE  RELATIVE  MACHINE  PRECISION.  IF  F  IS  A  UNIMODAL 
C  FUNCTION  AND  THE  COMPUTED  VALUES  OF  F  ARE  ALWAYS  UNIMODAL  WHEN 
C  SEPARATED  BY  AT  LEAST  EPS*ABf(X)  ♦  (TOL/3),  THEN  FHIN  APPROXIMATES 
C  THE  ABCISSA  OF  THE  GLOBAL  MINIMUM  OF  F  ON  THE  INTERVAL  AI,BX  WITH 
C  AN  ERROR  LESS  THAN  3*EPS*ABS  (FHIN)  +  TOL.  IF  F  IS  NOT  UNIMODAL, 
C  THEN  FMIN  MAY  APPROXIMATE  A  LOCAL,  BUT  PERHAPS  NON-GLOBAL,  MINIMUM  TO 
C  THE  SAME  ACCURACY. 

C  THIS  FUNCTION  SUBPROGRAM  IS  A  SLIGHTLY  MODIFIED  VERSION  OF  THE 

C  ALGOL  60  PROCEDURE  LOCALMIN  GIVEN  IN  RICHARD  BRENT,  ALGORITHMS  FOR 
C  MINIMIZATION  WITHOUT  DERIVATIVES,  PRENTICE  -  HALL,  INC.  (1973). 

C 

DOUBLE  PRECISION  A,B,C,D,E, EPS ,XM , P, Q, R, TOL1 , TOL2, U,V, W 
DOUBLE  PRECISION  FU,FV,FW,FX,  X 
DOUBLE  PRECISION  DABS ,DSQRT ,DSIGN 

C  IS  THE  SQUARED  INVERSE  OF  THE  GOLDEN  RATIO 

C  *  0.  5D0*(3.  -  DSQRT  (5.OD0)  } 

C 


EPS  IS  APPROXIMATELY  THE  SQOAHE  ROOT  OP  THE  RELATIVE  MACHINE 
PRECISION. 


EPS  =  1.ODO0 
10  EPS  =  EPS/2. 0D00 
TOL1  =  1.0D0  *  EPS 
IF  (TOL1  .GT.  1.0D00)  GO  TO  10 
EPS  =  DSQET  (EPS) 

:  INITIALIZATION 

A  =  AX 
B  =  3X 

*V  »  A  ♦  C*  (B  -  A) 

W  =  V 
I  =  V 
E  =  0.0  DO 
FX  =  F  (X) 

FV  =  FX 
FW  =  PX 

:  MAIN  LOOP  STARTS  HERE 

20  XM  =  0.5D0*  (A  +  B) 

TOL1  »  EPS*DABS (X)  ♦  TOL/3.0D0 
TOE 2  *  2. ODO*TOL 1 

:  CHECK  STOPPING  CRITERION 

IF  (DABS  (X  -  XM)  .LE.  (TOL2  -  0.5D0*(B  -A)))  GO  TO  90 

:  IS  GOLDEN- SECTION  NECESSARY 

IF  (DABS  (E)  .LE.  TOL1)  GO  TO  40 

:  FIT  PARABOLA 

P  *  (X  -  W)  *(PX  -  FV) 

Q  =  (X  -  V)  *<FX  -  FW) 

P  =  (X  -  V)  *Q  -  (X  -  W)  *R 
Q  =  2.  ODOO*  (Q  -  R) 

IF  (Q  .GT.  O.ODO)  P  =  -P 
Q  *  DABS (Q) 

R  *  E 

E  *  D 

:  IS  PABABOLA  ACCEPTABLE 

30  IF  (DABS  (P)  -GE.  DABS  (0. 5D0*Q*R) )  GO  TO  40 
IF  (P  .LE.  Q*  ( A  -  X))  GO  TO  40 


IP  (P  .GE.  Q*  (B  -  X))  GO  TO  40 

A  PAPABOLIC  INTERPOLATION  STEP 

D  =  P/Q 
U  =  X  +  D 

F  MOST  NOT  BE  EVALUATED  TOO  CLOSE  TO  AX  OB  BX 

IF  (  fO  -  A)  .LT.  TOL2)  D  =  DSIGN  (TOL1  ,  XH  -  X) 

IF  (  (B  -  0)  .LT.  TCL2)  D  =  DSIGN(TOL1,  XS  -  X) 

GO  TO  50 

A  GOLDEN-SECTION  STEP 

40  IF  (X  .  GE.  XM)  E  -  A  -  X 
IF  {X  .LT.  XU)  E  -  B  -  X 
D  =  C*E 

F  BUST  NOT  BE  EVALUATED  TOO  CLOSE  TO  X 

50  IF  (DABS  (D)  .GE.  TOL1)  0  =  X  ♦  D 

IF  (DABS  (D)  .LT.  TOL1)  0  *  X  +  DSIGN  (TOL1,  D) 

PO  =  F  (D) 

UPDATE  A,  B,  V,  W,  AND  X 

IF  (FO  .GT.  FX)  GO  TO  60 
IF  (0  .GE.  X)  A  »  X 

IF  (0  .LT.  X)  B  =  X 

V  =  W 
FV  =  FW 
W  =  X 
PW  *  PX 
X  *  0 
PX  =  FO 
GO  TO  20 

60  IF  (0  .LT.  X)  A  =  0 

IF  (0  .GE.  X)  B  *  0 

IF  (FO  .LE.  FW)  GO  TO  70 
IF  (W  .EQ.  X)  GO  TO  70 
IF  (FO  .LE.  FV)  GO  TO  80 
IF  (V  .EQ.  X)  GO  TO  80 
IF  (V  .EQ.  W)  GO  TO  80 
GO  TO  20 
70  V  *  W 
FV  *  FW. 

?  »  0 
PW  *  FO 
GO  TO  20 


END  OF  SAIN  LOOP 


90  FSIN  =  X 
RETORN 
END 


********************* ******************  *****************************  ** 

PROGRAM  NAME  :  MAPGRAF 

WRITTEN  BY  :  J. F.  BRANDEAU 

COMPILER  (S)  :  WATFIV  (DOUBLE  PRECISION) 


PURPOSE  :  USE  COEFFICIENTS  X  TO  PRODUCE  DATA  POINTS  THAT  WILL  BE  READ 
BY  SAS  TO  PRODUCE  A  3-D  PLOT  OF  A  STRESS-STRAIN  SURFACE.  PROGRAM 
CREATES  A  GRID  OF  STRAIN  VS.  LOG  STRAIN  RATE,  AND  USES  SUBROUTINE 
ZEROIN  TO  SOLVE  FOR  THE  VALUE  OF  STRESS  FOR  EACH  UNIQUE  COMBINATION. 
THE  DENSITY  OF  THE  GRID,  THE  RANGE  OF  THE  GRID,  AND  THE  NON- 
PERMISSIBLE  REGION  ARE  ALL  CONTROLLABLE.  THE  NON -PERMISSIBLE  REGION 
OF  THE  GRID  WILL  HAVE  NO  POINTS  ON  IT,  TO  PREVENT  EXTRAPOLATION  PAST 
EXPERIMENTAL  LIMITS.  THE  SHAPE  OF  THE  EDGE  IS  CALCULATED  USING  A 
SPLINE  CURVE  FIT  TO  THE  MAX  VALUES  OF  STRAIN  FOR  EACH  OBSERVED  STRESS 
THIS  IS  DONE  BY  SUBROUTINE  SPLINE  AND  SEVAI. 

VARIABLES  : 

A1  E  B1  :  A 1  IS  THE  MAX  OBSERVED  VALUES  OF  STRAIN,  ONE  FOR  EACH  VALUE 
OF  STRAIN  RATE  B1.  B1  MUST  START  WITH  LOWEST  VALUE  OF  STRAIN  RATE 
FIRST,  OP  THE  SPLINE  WILL  BE  INCORRECTLY  CALCULATED. 

K  t  THE  NUMBER  OF  UNIQUE  STRAIN  PATES  OBSERVED,  ALSO  THE  NUMBER  OF 
ELEMENTS  I  EACH  OF  A1  5  Bl. 

SIGMAX  :  GREATER  THAN  THE  EXPECTED  MAX  VALUE  OF  STRESS  (ESI)  TO  BE 
USED  AS  THE  UPPER  LIMIT  OF  SEARCH  FOR  ZEROIN.  FOR  EACH  UNIQUE 
STRAIN  RATE  THIS  MAY  BE  CORRECTED  DOWNWARD  BY  THE  PROGRAM  TO  PREVENT 
OVERFLOWS.  THIS  IS  SET  BY  THE  USER. 

A,B,N,D,C  6  COEFF  :  A ,B,N,D,  6  C  ARE  THE  MANTISSA  OF  THE  COEFFICIENTS 
FOR  THE  EQUATION,  AND  COEFF  IS  THE  EXPONENTS.  THE  PRODUCTS 
A  *  COEFF  (1),  B  *  COEFF (2)  ,  ETC...  ABE  THE  COEFFICIENTS. 

XI  :  THE  STARTING  POINT  FOR  VALUES  OF  STRAIN.  THIS  MUST  BE  GREATER 
THAN  ZERO,  AS  THE  EQUATION  IS  INDETERMINATE  9  STRAIN  =  0.0. 

Y1  :  THE  STARTING  POINT  FOR  LOG  RATE.  Y1  »  -3  IS  A  STARTING  POINT  OF 
STRAIN  RATE  =  0.001  /  SEC. 

DY  :  STEP  LENGTH  IN  LOG  RATE  (OR  Y)  DIRECTION. 


DX  :  * 


"  STRAIN  (OR  X) 


YHAX  :  MAX  VALUE  OF  LOG  RATE  TO  BE  USED  &  IS  ASSIGNED  BY  THE  USES. 

THOLD  *  MAX  VALUE  OP  STRAIN  TO  USE,  5  IS  SET  BY  PROGRAM  TO  EQUAL 
SAX  VALUE  IN  A 1. 


*00000010 

*00000020 

*00000030 

*00000040 

*00000050 

*00000060 

*00000070 

*00000080 

*00000090 

*00000100 

*00000110 

*00000120 

*00000130 

*00000140 

*00000150 

*00000160 

*00000170 

*00000180 

*00000190 

*00000200 

*00000210 

*00000220 

*00000230 

*00000240 

*00000250 

*00000260 

*00000270 

*00000280 

*00000290 

*00000300 

*00000310 

*00000320 

*00000330 

*00000340 

*00000350 

*00000360 

*00000370 

*00000380 

*00000390 

*00000400 

*00000410 

*00000420 

*00000430 

*00000440 

*00000450 

*00000460 

*00000470 

*00000480 

*00000490 

*00000500 


TOL  :  CONVERGENCE  CRITERIA  FOR  ZEROIN 


IX  :  t  OP  POINTS  ON  STRAIN  (X)  AXIS. 

IY  :  "  "  "  "  LOG  RATE  (Y)  " 

THE  TOTAL  GRID  WILL  CONSIST  OP  IY  *  IX  -  NON-PERHISSIBLE  POINTS. 

C1,D1,E1  :  WORK  VECTORS  FOR  SPLINE  5  SSVAL.  BUST  NOT  BE  CHANGED. 

I/O  REQUIREMENTS  : 

FILE  #5  :  MANTISSA  OF  COEFFICIENT  VECTOR  AND  EXPONENT  OF  COEFFICIENT 
VECTOR  ON  TWO  RECORDS. 


OPTIONS  :  NONE 


****************  *****41***4141***  *********************************  ******* 

IMPLICIT  REAL  *  8  (A-H,  N,  O-Z) 

EXTERNAL  FUNCT 

DIMENSION  Y  (40)  ,  X (40)  ,  YE  (40)  ,  COEFF(S) 

DIMENSION  A'l  (6)  ,  B 1  f *8)  ,  Cl  (6) ,  Dl(6),  El  (6) 

COMMON  X,  Y,  I,  J,  A,  B,  C,  D,  N,  8 ATEB,  RATED,  TLOG 
COMMON  /SUB  1/  GEPS 

DATA  A1/6.90D— 3,  6.5D-3,  5.7D-3,  5. ID-3,  10.5D-3,  9.2D-3/ 

DATA  B 1  /.001,  0.10,  10.0,  150.0,  0.0,  0.0/ 

SET  PROGRAM  PARAMETERS 

K  *  4 

YMAX  =2.3 
SIGMAX  =  25.0D0 
XI  -  0.05D-3  ;  Y1  *  -3.0D0 

TOL  =  1.0D-8  ;  IX  =  30  ;  IY  =  30 

THOLD  =0.0 
DO  5  J  =  1,  K 

IF  (A1  (J)  .GT.  THOLD)  THOLD  =  A1(J) 

5  Bt  (J)  =  DLOGIO  (B 1  (J) ) 


CALL  SPLINE  (K,  B1 
READ  (5,*)  A,  B,  N 
READ  (5,*)  COEFF 
A  =  A  *  COEFF  ( 1) 

D  =  D  *  COEFF  (4) 
PRINT,  »A  5 

PRINT,  *D  * ,D  ; 


A 1 ,  Cl,  D1,  El) 

D,  C 

B  =  B*  COEFF (2) 
C  =  C  *  COEFF  (5) 
PRINT,  *B  * ,  B 
PRINT , * C  »,C 


N  =  N  *  COEFF  (3) 
PRINT,  * N  * ,N 


CHECK  FOR  VALUES  THAT  WILL  CAUSE  OYER/UNDER  FLOW 


IF  (A  .EQ.  0.0D0)  THEN  DO 
TLOG  *  -80.0D0 


*00000510 
*00000520 
*00000530 
*00000540 
*0000055 0 
*00000560 
*00000570 
*00000580 
*00000590 
*00000600 
*00000610 
*00000620 
*00000630 
*00000640 
*00000650 
*00000660 
*00000670 
00000680 
00000690 
00000700 
00000710 
00000720 
00000730 
00000740 
00000750 
00000760 
00000770 
;0000780 
00000790 
00000800 
00000810 
00000820 
00000830 
00000840 
00000850 
00000860 
00000870 
00000880 
00000890 
00000900 
00000910 
00000920 
00000930 
00000940 
00000950 
00000960 
00000970 
00000980 
00000990 
00001000 


ELSE  DO 

TLOG  =  DIOGIO  (A) 

ENDIF 

DX  *  (THOLD  -  XI)  /  DFLOAT  (IX-1) 

DT  *  (YHAX  -  71)  /  DFLOAT  (IY-1) 

DO  10  J  =  1,  17 

'  Y(J)  *  DFLOAT  (J-1)  *  DY  ♦  Y1 
10  YE  (J)  =  10.0  **  Y(J) 

DO  20  I  =  1,  IX 

20  X(I)  =  DFLOAT  (1-1)  *  DX  ♦  XI 

CALCULATE  MACHINE  EPSILON  FOE  ZEROIN 

GEPS  *  1.0D0 

24  GEPS  =  GEPS/2.0D0 
TOL1  »  1.0D0  ♦  GEPS 

IF  (TOL1  .GT.  I.ODO)  GO  TO  24 

BEGIN  SOLUTIONS 
DO  40  J  =  1,  IY 

CORRECT  BX  IF  NEEDED  TO  PREVENT  OVERFLOW 

\ 

BI  a  SIGMAX 

RATEB  =  DLOGIO  (YE  (J)  **  B) 

25  AX  =  N  *  DLOGIO  (BX)  +  RATEB 

IF  (AX  .GT.  75.0D0)  THEN  DO 
BX  a  BX  -  0.5D0 
GO  TO  25 
ENDIF 

RATEB  *  YE  (J)  **  B 
RATED  =  C  *  YE ( J)  **  D 
Z  *  X ( 1 )  *  RATED  /  1.3 

TEMP  =  SEVAL  (K,  Y  ( J)  ,  BI,  A1,  Cl,  D1,  El) 

ALLOW  FOR  SPLINE  RIPPLE 

IF  (TEMP  .GT.  (THOLD  *  0.97D0)  )  TEMP  *  THOLD 

DO  30  I  a  1,  ix 

IF  (T(I)  .GT.  TEMP)  GO  TO  40 
AX  =  Z 

Z  «  ZEROIN  (AX,  BX,  FUNCT,  TOL) 

WRITE  (6,400)  X(I),  Y(J)  ,  YE(J),  Z 
30  CONTINUE 
40  CONTINUE 
STOP 

400  FORMAT  (1H  ,5(1PD13.5)) 


00001010 
00001020 
00001030 
00001040 
00001050 
00001060 
00001070 
00001080 
00001090 
00001100 
00001110 
00001120 
00001130 
00001140 
00001150 
00001160 
00001170 
00001180 
00001190 
00001200 
00001210 
00001220 
00001230 
00001240 
00001250 
00001260 
00001270 
00001280 
00001290 
00  001300 
00001310 
00001320 
00001330 
0^01340 
OOuJ1350 
00001360 
00001370 
00001380 
00001390 
00001400 
00001410 
00001420 
00001430 
00001440 
00001450 
00001460 
00001470 
00001480 
00001490 
00001500 


END 

DOUBLE  PRECISION  FUNCTION  FUNCT (STRESS) 

IMPLICIT  REAL  *  8  (A-H,  N,  O-Z) 

DIMENSION  Y  (40)  ,  X{40) 

COMMON  X,  Y,  I,  J,  i,  B,  C,  D,  N,  RATEB,  RATED,  TLOG 
T1  »  STRESS  /  RATED 
IF  (STRESS  .GT.  1.0D-1)  GO  TO  10 
IF  (  (N  *  DLOGIO  (STRESS) )  .GT.  -60.0DO)  GO  TO  10 
HOLD  =  -25.0DO 
GO  TO  15 

10  T2  =  STRESS  **  N 
HOLD  =  DLOGIO  (T2) 

15  IF  ((HOLD  ♦  TLOG)  .LT.  -17.0D0)  THEN  DO 
■  FUNCT  =  T1  -  X(I) 

ELSE  DO 

FUNCT  =  T1  ♦  A  *  T2  *  RATEB  -  X(I) 

ENDIF 
25  RETURN 
END 

SOPTIONS  NOLIST 

DOUBLE  PRECISION  FUNCTION  ZEROIN (AX,BX,FrTOL) 

DOUBLE  PRECISION  AX,BX,F,TOL 

DOUBLE  PRECISION  A, B, C,D , E, EPS ,FA, FB , FC.TOL1 , XM ,P ,Q ,R ,  S 
DOUBLE  PRECISION  DABS,DSIGN 
COMMON  /SUB  1/  EPS 
A  -  AX 
B  *  BX 
FA  =  FfA) 

FB  =  F  ( B) 

20  C  ■  A 
FC  =  FA 
D  =  B  -  A 
E  ■  D 

30  IF  (DABS  (FC)  .GE.  DABS  (FB)  )  GO  TO  40 
A  =  B 
B  =  C 
C  =  A 
FA  *  FB 
FB  *  FC 
FC  *  FA 

40  TOL1  =  2. 0D0*EPS*D ABS (B)  ♦  0.5D0*TOL 
XM  '»  .  5*  (C  -  B) 

IF  (DABS  (XM)  .LE.  TOL1)  GO  TO  90 
IF  (FB  .EQ.  0.0D0)  GO  TO  90 
IF  (DABS  (E)  .LT.  TOL1)  GO  TO  70 
IF  (DABS (FA)  .LB.  DABS  (FB) )  GO  TO  70 
IF  (A  .NE.  C)  GO  TO  50 
S  *  FB/FA 
P  *  2.0D0*XH*S 
Q  *  1.0D0  -  S 


00001510 

00001520 

00001530 

00001540 

00001550 

00001560 

00001570 

00001580 

00001590 

00001600 

00001610 

00001620 

00001630 

00001640 

00001650 

00001660 

00001670 

00001680 

00001690 

00001700 

00001710 

00001720 

00001730 

00001740 

00001750 

00001760 

00001770 

00001780 

00001790 

00001800 

00001810 

00001820 

00001830 

00001840 

00001850 

00001860 

00001870 

00001880 

00001890 

00001900 

00001910 

00001920 

00001930 

00001940 

00001950 

00001960 

00001970 

00001980 

00001990 

00002000 


GO  TO  60 
50  Q  =  FA/FC 
B  »  FB/FC 
S  =  FB/FA 

P  =  S*  (2.0D0*XH*Q*  (Q  -  E)  -  (3  -  A)*(H  -  1.0D0)) 

Q  =  (Q  -  1 . 0D0)  *  (H  -  1.0D0)  *  (S  -  1.0D0) 

60  IF  (P  .GT.  0.0D0)  Q  =  -Q 
P  =  DABS  (P) 

IF  (  (2  . 0D0*P)  .GE.-  (3.0D0*XN*Q  -  DABS  (TOL 1*Q)  ) )  GO  TO  70 
IF  (P  .  GE.  DABS (0. 5D0*E*Q) )  GO  TO  70 
E  =  D 
D  =  P/Q 
GO  TO  80 
70 • D  *  XH 
E  =  D 
80  A  =  B 
FA  *  FB 

IF  (DABS  (D)  .GT.  TOL1)  B  *  B  *  D 

IF  (DABS  (D)  .LE.  TOL1)  8*8  +  DSIGN(T0L1,  IB) 

FB  *  F (B) 

IF  (  (FB*  (FC/DABS  (FC)  )  )  .GT.  O.ODO)  GO  TO  20 
GO  TO  30 
90  ZFROIN  *  B 
RETURN 
END 

SUBROUTINE  SPLINE  (N,  X,  Y,  B,  C  D) 

INTEGER  N 

DOUBLE  PRECISION  X  (N)  ,  T  (N)  ,  B  (N)  ,  C(N),  D  (N) 

INTEGER  NM1,  IB,  I 
DOUBLE  PRECISION  T 
HB1  =  N-1 

IF  (  N  -LT.  2  )  RETURN 
IF  (  N  .LT.  3  )  GO  TO  50 
D  (1)  *  X  (2)  -  K1) 

C  <21  *  (1(2)  -  Y(1))/D(1) 

.  DO  10  I  =  2,  N51 

d  (i)  »  x(i+n  -  x(D 

B  (I)  *  2.*{D(I-1)  ♦  D  (I) ) 

C  (1  +  1)  =  (Y  (1+1 )  -  Y(I))/D(I) 

C  (II  *  C  (I ♦  1 )  -  C  (I) 

10  CONTINUE 

B  (1)  a  — D  (  1  ) 

B  (N)  a  -D(N-1) 

C(1)  *  0. 

C (N)  *  0. 

IF  (  N  .EQ.  3  )  GO  TO  15 

C  (1)  *  C  (3)  /(X  (4)  -X  (2) )  -  C  (2)  /  (X  { 3)  -X  (1) ) 

C  (N)  *  C(N-1)/(X(N) -X(N-2)  )  -  C  (N-2)  /  (X  (N-1) -X  (N-3)  ) 

C(1)  *  C  ( 1)  *D  { 1)  **2/  (X  (4)  —X  (1) ) 

C (N)  »  -C(N)*D(N-1)**2/(X(N)-X(N-3)) 


00002010 

00002020 

00002030 

00002040 

00002050 

00002060 

00002070 

00002080 

00002090 

00002100 

00002110 

00002120 

00002130 

00002140 

00002150 

00002160 

00002170 

00002180 

00002190 

00002200 

00002210 

00002220 

00002230 

00002240 

00002250 

00002260 

00002270 

00002280 

00002290 

00002300 

00002310 

00002320 

00002330 

00002340 

00002350 

00002360 

00002370 

00002380 

00002390 

00002400 

00002410 

00002420 

.00002430 

00002440 

00002450 

00002460 

00002470 

00002480 

00002490 

00002500 


15 

DO  20  I  =  2,  N 

'1 

00002510 

T  =  D  (I-1J/B  (1-1) 

00002520 

B  (I)  =  B  (I)  -  T*D  (1-1) 

00002530 

C  (I)  =  C  (I)  -  T*C  (1-1) 

00002540 

20 

CONTINUE 

00002550 

C  (N)  =  C(N)/8(N) 

00002560 

DO  30  IB  =  1,  NH1 

00002570 

I  =  N-IB 

00002580 

C(I)  =  (C(I)  ^  D{I)*C(I  +  1))/B(I) 

00002590 

[■  30 

CONTINUE 

00002600 

B  (N)  =  (Y(N)  -  I  (Nf!1)}/D(NH1)  ♦  D  (  NH1 )  *  (C  (N(1 1 )  ♦  2.*C[N)) 

00002610 

DO  40  I  =  1,  NH1 

00002620 

■ 

B  (I)  =  ( Y  (1  +  1)  -  Y(I))/D(I)  -  D  (I)  *  (C  (I+1)  +  2.  *C  (I)  ) 

00002630 

D  (I)  =  (C(I+1)  -  C(T))/D(I) 

00002640 

C(I)  =  3.*C(I) 

00002650 

40 

CONTINUE 

00 002660 

C  (N)  =  3-  *C  (N) 

00002670 

D  (N)  =  D(N-1) 

00002680 

RETURN  ^ 

00002690 

50 

B  (1)  =  (Y(2)-Y(1))/(X(2)-X(1)) 

00002700 

c  ( 11  =  0.  • 

00002710 

D  { 1 )  =  0. 

00002720 

B  (2)  =  B(1) 

00002730 

C (2)  =  0. 

00002740 

D(2)  =  0. 

00002750 

RETURN 

00002760 

END 

00002770 

DOUBLE  PRECISION  FUNCTION  SEVAL(N,  Ur  X,  Y,  B,  C,  D) 

00002780 

INTEGER  N 

00002790 

DOUBLE  PRECISION  U,  X  (N)  ,  Y  (N)  ,  B  (N)  ,  C  (N)  ,  D  (N) 

00002800 

INTEGER  I,  J,  K 

00002810 

DOUBLE  PRECISION  DX 

00002820 

DATA  1/1/ 

00002830 

IF  (  I  .GE.  N  )  I  *  1 

00002840 

IF  (  U  -LT.  X(I)  )  GO  TO  10 

00002850 

IF  (  U  .LE.  X(I*1)  )  GO  TO  30 

00002860 

10 

I  =  1 

00002870 

J  *  'N+  1 

00002880 

20 

K  =  (I  ♦  J)  /2 

00002890 

IF  {  U  -LT.  X(K)  )  J  ■  K  - 

00002900 

IF  (  U  -GEi  X(K)  )  I  =  K 

00002910 

IF  (  J  -GT.  1  +  1  )  GO  TO  20 

00002920 

30 

DX  »  0  -  X  (I)  v  ...  . 

00002930 

SEYAL  =  1(1)  ♦  DX*  (B  (I)  ♦  DX*  (C  (I)  ♦  DX*D(I))) 

00002940 

RETURN 

00002950 

j 

END 

00002960 

rr^rr 


:  soft  i  on  — - - — . .  • — - “ — - — - 

.  PPOGPAM  NAME:  STRESS  . .  ~  — 

REAL  DUMMYC24)  ,  CIUTCM,  XC15),  L95 (4) ,  YO I ( 15 T .  '  ~  ~  ' 

PEAL  a,  B,  C,  CG(  15)  ,  WETGHTfrsrv  rN6RTf45),  DIT4) 

REAL  DM  EGA  (  3)  »  FJ  (  3 ) ,  T  JT3  )  ,  ALPHA  (31  ,  ACCEL C 3) ,  U2f3 > " 

T  REAL  POSIT(74T;.~CWSTKNT^r” ’ — — r - - - - *~T • — ~ — - 

■  r  NTE  3  c  h  S  3  J  7th  E CRT 47, M njNT“lT0^  1  V FT ,"RT57~3T~^'i  TETST^TT 

INTEGER  OUO  - - - —  ~ 

'  •  loilivaijljj^  ^FORCEfl  1  ^  (  S  2,'  TEMP  ,  SUBJH _ _ 

COMMON  / 06/  B0NEPC8),  80NEA18),  BONE! 181,  LlL2r  WIH2,  Al,  I NEP l > 

5  '  R0RI3 ,  TNEP2,  RO"  '  '  -  '  '  . . 

--WHS 

DATA  GINCH  /  386.0836  /,  PI  /  3.14T593  /,  W95  /  21T.0  / 

'  DET  f  A  ,  a  1  =  SQP.T  f  A  *  A  •*■  B  *  B  /  2  ~T~ . .  ‘  * 

WAX  »  0.0  - - -  -  -  . . • -  -  — .  - 

‘  WRITE  l  3,  5) 


000000  10 
00000020  _ 
00000030 

:SM§F- 

00000060 

00000070 


nroo‘oowo 

00000100 
000001 10 

mmF- 

OOOOOL4C 

00000150 

mm= 

00000180 

00000190 

00000220 
00000230 
“05JODO  240LZ 


WRITE  (8,61  1,  SEGTI) -  - - -  - - - - -  - 

~  6*  FORMAT  (IH-  ,  13,  1  r*  A 4 V"  “ - - - - - - - - - - - 

”  T-  r-n WT I  MTTF ~ •  '^u  '  "  „:m„.  I.,,... ..  .  .....  ■  -a-»— — — — 

- ^ ^  •«>»..  .till..  .m  Mint  i.  i.  |»  .  tt  .» . .,  ...  ,,  i 

~  PE  ADS (  5 1'^MTV  DLTt~  'KCON4'r  '  1° JSIT( TlT'^r  =“~T, ~KC0N4 1 - — * 

— -  — YpC  ^  *rr  ~~  ”'wp1"1  *'*'  w  <— u W |>i  'WT'  WiyWTK'i  M  M 

- D J  ‘3~T  l,",< — 

Til  IF  IPOSlTCn  VNE.  J>  GO  TCT'8;~ ~~  — . •T;~irp-" _ -"T, - -r  - - 

l  XJ N0~~MASS,  D1~ANQ  ~StM  1- MAJOR.'  A  XES~~0F~~ SEGMENT  J~ FROM  FILE  RS  ~~ - 

'  — r irApr^i r fr~'‘  — ;-r=-  a  ■.■  -- 

IF  CJ  ,NE.  II  TALL  FOP  TJ-1,  53 - —  - -  - 

READ'  (5)  <»_  MASS,  A  ,  8,  C  V  -  “  '  - 

-  f  huN’o~  \T1 

WRTTE'  (  3,  150)  ‘  J  ,  “SE3('J)T'NT T'NSEGS' - -  -  — ' 

zzzzmz  t-  zsissssa 

' . _J.aii.JfES.  *  3  JiT.-T  "'T'T - ~~7Z . - I - j  - - 


00000260 

00000270 

WHF1 

00000300 

00000310 


“00000340 

00000 35D 

zgagiuigeTr- 

00000  3  7  0 
-00000380 
00000 3°0 

raaaffl34Qoi_ 
00000410 
0000042  0 
000004’ : 
'JOQSDQS^ll. 
00000450 
.00000460 
00000470 
3K30CEK20Z: 

000004-0 

00000500  _ 


✓  AD-AU2  AM  DUKC  UN IV  DURHAM  NC  -DEPT  Of  MECHANICAL  EN6XNECRIN6  A— ETC  F/6  6/19  \ 

DEVELOPMENT  OP  A  STRAIN  RATE  DEPENOEMT  L0N6  BONE  INJURY  CRXTERI— ETC<U> 

JAN  62  T  K  HI6MT  APOSR-01-0062 

UNCLASSIFIED  AFOSR-TR-62-01M  NL 


non  non  :  ;  t  r>no< 


IF  (YES  .r,T.  0)  WRITS  ,(  3,  156)  (  PQS.TTtJ  >,  !  »  IL,  IRI 
S ST  JNin  TO  tluAL  PROPER  COUNTER 'I N  GEOMETRY  FILE'  45  “ 


'  R  5aO^S7*1~Ch  tCK 

go  to  no,  20,  ja,  2ff,;3ar  4o*:3o,:4orii:j: 

10  JNUM  =1 

rrn **T"  upwwirv-T'^ti  ******** '  -ty  n  5 v« 

- ug J21 *  - 1  ,  .  ,JaLS3«>£. £»5S5S  iBS 

GO  TO  50 

20  JNUM  =2  •  - - '"r  •  - 

CD”)  =  -1  .  .  -  -  ~~~~ 

. ..  ~  calc,  fop  (CHEttcngnT^Tpygass^? 

GO  To  *50 

30  JNUH  -  3  . ZZZZZZIZZZ  - .~-:"TT — 

ODD  I  - - 


trilfl ■  « m  mJijAwii  3» 


GO  TO  50 

40  JNUW  =  4 
~  ODD  =  -I 

r v 


viiflus  y 


1  =  A  l 

X2  *  XI  -*  XI 
X?  =  X2  *  XI 


00000510 

•00000520 

00000530 


.00000560 

00000570 

□joflcossa: 

00000590 

00000600 

OOOOOfclO 

3o.aofl.62a. 

00000630 

00000640 

00000650 


(  jJlIUtT#  /VfJltl 


00000670 

00000680 

00000690 


■IhiM'IIH 


OOOOOTZO 
00000730 
00000740 
00000750“ 
00000760 
00000770 
OOTBO 
OOOO0T9O 


00000800 

_ 00000810 

r^rrTTooaoofzo"* 

0TH30 


00000880 

00000890 


mi  1 1 1«!^ —— — 

. >■  <? .**+»*. n— kOl*  lu*  *l>i< 

PI  III,  — — — — — — — — — —  ilililili I  —j 


00000910 

00000920 

00000930 

mmz 

00000960 

00000970 


F IKO  HEtGHT  fOR  WASS r  OF~:F8FE— BOOr’SECTrOfT 


non  nnnrtn  o  o<~>o  o<~>r> 


60  CONTINUE 


AVAILABLE  (CHECMJNUM) 

HOLLOW  TUBE  FOR  .BONE  T CHECK f 4 NUN I 


CONSTANT  CROSS-SECTION. 


1 2  -«  49  I  m  aril*’  di  9 


00001180 
0000 1190 

■  0O00L2.00. 


:poss-sect:on  00001220 

_ _  00001230 

■ - - - 00001240 


00001260 

00001270 

00001280 


[%toats 


00001520 
1530 


.540 


STOP 


95  OCT  .145  KK  *  T*  NT 


0000.15  50 


00001600 


c  .:f:ino  a 


READ  C5I  TIME 

IF  (J  .  N  E  .  1 )  CAt L  TO fi  l  J-1 


*  r*fi: 


00001680 

00 

0000 


I'UhIL 


0  L 

0A 

OS' , 6  F 12 

TO" 

96 

HEP  CON 


C  BEGIN  LOOP  FOR  EACH  FREE— BODY  SECTION  HT  -  - 

C  ..  C ALCULATE  STKE5SE S  ON  EACH  FREE-BODY  SECTI  ON  US ING  MOHR1  S~~C  IRC L 

***  •  rS?f  *  W  .  •*  .'  "IMP 


C— AMCLSXATIC  ;EgfILI - 


100  WHITE  ( 3, Z50)  . . ' 

_ -0Q145  i  =  i;.:nsegs 
.  BTrjr^ x^tu  ifr\yryr^r~ 


00001800 

00001810 

00001820 


II I » W.fi. 


mm 


p*9*cn&* 


1'ZZ^Z.lS.  -Lfroooo 

00 


00  105  K. *_.! 


C  F I  NIT  JIOMENTS-  ANP-^UM  „.FOFCES-DUE;XO  _’»aTHEF.:CONSTRATNr*FORCES":  .... 


00001920 

00001930 


MUM— II 'i^Mlih  Will  Will 


IF  ( XPu I  NT  i  107,  107,  106 
106  MOMENT  (l)j».  POSIT  (YES  +2  1  ♦  CNSTRNI3? 
M0MENTJ.2)  *  XPOINT  *  ‘ 


»!., 4  <**!:<  WT^OY-j  •  >«.  » -if!  -*r>«  Jit*  s»  4 

"  •Ii 


**  «  TNT  7*4 


n^*~***  *v**>^^  y.  .  y»-  k  ^ 4tc-r 


*44******* >*■  *5«»xw.  -jA. •«..? 


>IJO 


C  SOLVE  MOHS  • 

C 

DO  125 . K2 


OR  EACH  END  OF  DIAMETER  OF  TU9E 


T  5  mP^DcT  (TA  uv  Su 
TEMP2  *  SXGIK2)  f  2. 
SA  '  =  TEM02  TE  MP~ 


00002*560 
000025 
000025 


mm 


125  coin  .  _  _  _  _____ _ 

=  AHA XI  rocm  i  r,~0UTT4TT 


,*i»3»r4*^s.»4S  i-T 


SlMAXJLGO_.Tp_iJ0. 


50002640 
■'"■  0000265CT 
. .  0000  2660 


HwSilaMBtNffMi  *  nil !  »*U.T*U 


00002680 

Q0QQ2690 

00002700 

rz^mqgpzTia; 

00002720 

.  "00002730 

- 00002740 


XT  »~XTII  .  ...  • 

140  wRTTg  r srrao roerx mirxrTTvnjirr 

C9  •  H«m  s~*  ” 

~r-;5"cn 

STOP 

C  . .  . 

ISOrFOTMA 


*  •  mc^t  AKc*»i3t*  nut 

*  ll»  *  FR  EE— BOOT  SECTIONS  MrtL  BF  USEOn-  r - ~ - - - 

15 5 “FOR MAT  (  FH-  f  T  SEGMENT  MASS' =  »,  FI0,3‘,  r  XB  r.  »‘/1'  '“SEMr-MAJO5_A'XES » 


V-nfinmmsmmmets  <wnm»x»  w  A^rfelu 

— - - - MiMsMaWil 


kUHUXUlUCtlVIllUMUMint 


S  •*  PROXIMAL  JOINT  TODISTAtJOfNTLENGTH  10V2, TMi  »7 

156  FORMAT  T r-tT OTHER  COWSTR  ATNT  FQRC5^~  LOCATTUN5^,lFro.'  3T~  ~  - 


imyilMPW  JLMBT  UltxlMtMHvi.OUU  *-*«* 

ICHEdLEl 


MMKlWWf' .  *  j  HUM*. 

MPJKlIMMMliaili'iauWBLgIM 


r  <rn»  >  *  «— >» 


[•HMtUIAl 


170  FORMAT  t  ♦-GEOHETRr  PATS" FOR  TMIS  LIMIT  — VARIABLE  CR0S5-SECTr0Nrr  00002970 
*  »-  HOLLOW  TURErA' -  BA  SED  ON  "EXP EKTH ENTTIT 0 ATS :  V - '  “  “  00002980 


AC  2 1  *;m23 C3LJ f 6 Ml 


fi*  Wfcjg:*  »»  #>.<»»,>  4-rv» 


00003220 

00003230 

00003240 


*wwM  !*'•  U^r4*tM 


